The burden of hand, foot, and mouth disease among children under different vaccination scenarios in China: a dynamic modelling study

Background Hand, foot, and mouth disease (HFMD) is a common illness in young children. A monovalent vaccine has been developed in China protecting against enterovirus-71, bivalent vaccines preventing HFMD caused by two viruses are under development. Objective To predict and compare the incidence of HFMD under different vaccination scenarios in China. Methods We developed a compartmental model to capture enterovirus transmission and the natural history of HFMD in children aged 0–5, and calibrated to reported cases in the same age-group from 2015 to 2018. We compared the following vaccination scenarios: different combinations of monovalent and bivalent vaccine; a program of constant vaccination to that of pulse vaccination prior to seasonal outbreaks. Results We estimate 1,982,819, 2,258,846, 1,948,522 and 2,398,566 cases from 2015 to 2018. Increased coverage of monovalent vaccine from 0 to 80% is predicted to decrease the cases by 797,262 (49.1%). Use of bivalent vaccine at an 80% coverage level would decrease the cases by 828,560. Use of a 2.0× pulse vaccination for the bivalent vaccine in addition to 80% coverage would reduce cases by over one million. The estimated R0 for HFMD in 2015–2018 was 1.08, 1.10, 1.35 and 1.17. Conclusions Our results point to the benefit of bivalent vaccine and using a pulse vaccination in specific months over routine vaccination. Other ways to control HFMD include isolation of patients in the early stage of dissemination, more frequent hand-washing and ventilation, and better treatment options for patients. Supplementary Information The online version contains supplementary material available at 10.1186/s12879-021-06157-w.


Introduction
Hand, foot, and mouth disease (HFMD) is a common infectious disease mainly caused by various enteroviruses. HFMD usually affects children under age of five, with a incidence rate of approximately 2400 cases per 100,000 in 2018 in China [1,2], and occurs more often in children under three [3]. After 2007, the HFMD epidemic was in an uptrend in China accompanied by serious outbreaks [4]. Although HFMD is usually self-limiting, it can result in complications associated with the central nervous system or death once progressing to severe cases [5,6]. There are more than 20 types of enterovirus leading to HFMD and Enterovirus 71 (EV71) and Coxsackie virus A16 (CV-A16) are the most commonly reported [3]. EV71 accounts for 70% severe HFMD cases and 90% HFMD-related deaths in mainland China [7].
Mathematical models of enterovirus transmission are increasingly being used to design and inform public health interventions to reduce the spread of HFMD. For example, Chuo et al. developed a deterministic susceptible-infectious-recovered (SIR) model to predict the number of infected people during an HFMD occurred in Sarawak, Malaysia and found that the disease spread quite rapidly and the parameter that may be able to control the number of susceptible persons [8]. Yang et al. used a SEIR model, adding a compartment to represent the exposed people (E) to capture the incubation period of HFMD and showed the spread trend of HFMD in China, and the basic reproduction number indicated that it was an outbreak [9]. By adding in a quarantined compartment (Q) in a SEIQR model, Liu simulated HFMD transmission and found that quarantine is the optimal strategy for HFMD protection [10].
EV71 vaccine is a new vaccine recently developed in China. The Vero cell-based EV71 inactivated vaccine has been shown to induce immune responses to EV71 among infants and young children between 6 and 59 months of age [11]. By March 2013, this inactivated vaccines had completed phase I-III clinical trials in Mainland China, showing good safety and greater than 95% protection against HFMD caused by EV71 [11][12][13]. Three manufactures in mainland China have recently completed phase III trials for EV71 vaccines and have been approved by China Food and Drug Administration in December 2015 [14]. In total, three monovalent EV71 vaccines have been licensed while bivalent EV71/CA16 vaccines are also under development and have been shown to induce protective immunity against both EV71 and CVA16 in vivo and in mice [15]. Further empirical studies are needed to assess the safety and effectiveness of these different vaccines in preventing HFMD.
As empiric data emerge on the safety and individual-level effectiveness of these vaccines, what remains unknown is their potential population-level impact on outbreaks and incidence of HFMD. Existing and future EV-71 and CV-A16 vaccines could differentially affect the incidence of HFMD. Thus, the aim of our study is to compare the predicted incidence of HFMD under different vaccination scenarios in China.

Methods
We developed and analyzed a deterministic, compartmental, mathematical model of enterovirus transmission to reproduce the observed HFMD epidemic in China from 2015 to 2018. The modelled population includes one age-group to reflect all children from age zero to 5 years of age.

Data sources Literature review
We conducted a comprehensive review of peer-reviewed and grey literature including government reports to obtain the following model parameters: transmission rate, duration of infectiousness, serotypes distribution, and vaccine efficacy.
National HFMD surveillance dataset HFMD refer to acute infectious diseases caused by enteroviruses such as human enterovirus 71 (EV-A71) and Coxsackie virus A group 16 (CV-A16), which are more common in preschool children. Severe cases are mostly caused by EV-A71 infection. The diagnostic criteria refer to the People's Republic of China Hygiene Industry Standard for Hand, Foot and Mouth Disease Diagnosis (WS 588-2018). We obtained the number of reported cases of HFMD per month from 2015 to 2018 in China from the official infectious disease report from Chinese Centre for Disease Control and Prevention (Table S1) [2]. Since January 1, 2008, all probable and laboratoryconfirmed HFMD cases have been reported to Chinese Centre for Disease Control and Prevention and by May 2, 2008, HFMD was a mandatory notifiable disease (III).

Demographic data
We obtained national-level estimates of the total number of children under 5, the birth rate, and the under-5 mortality rate from 2015 to 2018 from the National Bureau of Statistics and Provincial Bureau of Statistics [1].

Model description
We developed an extended SEIR 6/8-compartment model to evaluate the spread and incidence of HFMD. The population is divided into six or eight compartments as shown in Fig. 1, based on the natural history of HFMD and intervention with vaccination. Specifically, the total population (N) was divided into four categories: susceptible to HFMD (S), exposed to HFMD-related pathogens (E), infectious (I) and recovered (R). With the exception of the recovered compartment, the other three compartments (S, E, and I) were further divided into 2 or 3 compartments based on vaccination, type of pathogens and hospitalization: vaccinated (V); exposed (E s ); vaccinated and exposed (E v ); infectious but not hospitalized including mildly ill patients and asymptomatic (I N ); infectious hospitalized with pathogens other than EV71 or CVA16 (I 2 ); infectious hospitalized with EV71-HFMD or CVA16-HFMD (I 3 ) to capture the bivalent vaccine. Considering that the EV71 monovalent vaccine was actually introduced in 2017 in China, we should modify the meaning of I 2 and I 3 compartments in 2017 and 2018, where I −e h indicates people who are infectious with pathogens other than EV71 and I e h indicates people who are infectious with EV71-HFMD. The total population (N) is therefore given by N = S + V + E s + E v + I N + I 2 + I 3 + R. Fig. 1 and Table 1 describe the model parameters.
The dynamics of HFMD is described in the following set of coupled ordinary differential equations to represent the dynamics of HFMD without vaccination.
The following set of eqs.
(2) represent the dynamics of HFMD with the bivalent vaccine (protective against EV71 and CV-A16 simultaneously) and the monovalent vaccine (protective against EV71).
Formula (1) indicates the dynamic of HFMD under no vaccination. In formula (2), the vaccines both indicate the bivalent vaccine which can prevent children from being infected by EV71 and CV-A16 simultaneously and the monovalent vaccine which can only prevent children from being infected by EV71. We explicitly simulate the introduction of the EV71 vaccine in 2017 because the monovalent vaccine was introduced at that time in China (in December 2016) [14] .

Model assumptions
We assumed that: a) All the newborns are susceptible. b) The epidemic was established by 2015 with the initial number of infectious [I 2 (0) + I 3 (0)] on January 1, 2015 equal to the number of reported cases in January 2015 divided by 31. c) We assumed that a vaccinated person has received the entire vaccination procedure (two doses of vaccination) and will be protected against EV71 or CV-A16 immediately after completing the full course of inoculation with EV71 or CV-A16 vaccine, and there is no cross-protection effect with other virus subtypes. The possibility of dominant virus subtypes of HFMD changing after vaccination is not considered. d) There is no HFMD-attributable mortality. According to the HFMD report of Chinese Centers for Disease Control, the number of deaths due to HFMD is extremely small, about one in ten thousand [2]. e) Since the EV71 vaccines were launched in China in December 2016, the coverage of the vaccines is still very low. We assume that there is no impact of vaccination on the transmission rate β, and everyone was susceptible before and including 2016. f) We assumed contact patterns reflect mass action without heterogeneity by age nor mixing by agegroups. Thus, the contact rate is subsumed within the biological transmission probability β. g) We assumed I N include asymptomatic

Parameter estimation and model fitting
We estimate parameters by minimizing the sum of squares (MSS) [17] using the MATLAB R2018a (version 9.4) tool fminsearch [21] which is used to perform unconstrained nonlinear minimization. Posterior parameter values were obtained when the results of fminsearch converged.

MSS
Where I 2 indicates the number of no EV71 or CVA16 HFMD ðI −ec h ) / no EV71 HFMD ðI −e h Þ, I 3 indicates the number of EV71 or CVA16 HFMD ðI ec h Þ / only EV71-HFMD ðI e h Þ. In order to test the goodness of fit between our model and the observed data, we used a Chi-Square Test for the following hypotheses: Null hypothesis H 0 : The modeled results (I 2 + I 3 ) are equal to the observed number of HFMD cases shown in Table S1. Alternative hypothesis H 1 : The modeled results are not equal to the observed values shown in Table S1. The basic reproductive number (R 0 ) is of a great significance in Epidemiological studies. It is the expected number of secondary cases produced, in a completely susceptible population, after introducing a typical infective individual [22]. In this study, we can compute the basic reproduction number (R 0 ) in 2015-2018 using the method of van den Driessche and assuming that each year represents an independent epidemic with a fully susceptible population on January 1 of each year [23]. See the Additional file 2 for the derivation process and results. In this model, R 0 is related to transmission rate (β), progression rate (α), the proportion of I N , I 2 , I 3 (P 1 , P 2 , P 3 ), the recovery rate (γ 1 , γ 2 , γ 3 ), the mortality rate (μ) and the vaccine efficacy (VE).

Sensitivity analysis
Considering the formula of R 0 (Formula (2) in the Additional file 2, the expression of R 0 in 2015 and 2016), the basic reproductive number is affected by β, α, P 1 , γ 1 , γ 2 and γ 3 while α is fixed to 1. Therefore, we take a derivation with respect to β, P 1 , γ 1 , γ 2 and γ 3 , respectively. Derivation results are as follows: We further conducted a sensitivity analysis to measure the impact of these five parameters (β, P 1 , γ 1 , γ 2 and γ 3 ) on the number of EV71 or CV-A16 HFMD (I 3 ) cases using Vensim PLP software (version 7.3.5, Ventana Systems, Inc.)

Analyses of different vaccination type and coverage scenarios
We used the model to simulate the HFMD incidence of I 2 and I 3 under different vaccine coverage scenarios. To 210,400 person/day); In each of the above cases, the coverage rates of bivalent vaccine coverage were 0, 10, 30, 50, 80 and 100%. We estimated the rate of vaccination for a certain vaccination coverage level to be the number of children < 5 divided by 365 and multiplied by the coverage level. We assumed vaccine protection rate was 95% and there was no loss of protective immunity following vaccination.

Analyses of pulse vaccination
We then compared pulse vaccination strategies via increasing the vaccination rate in March, April and September (before the two peaks of HFMD incidence of a year) from 26,300 person/day (10% coverage rate), 78, 900 person/day (30% coverage rate), 131,500 person/day (50% coverage rate) and 210,400 person/day (80% coverage rate) to 1.2 times, 1.4 times, 1.6 times, 1.8 times and 2.0 times, respectively. The annual vaccine coverage rates under different pulse vaccination strategies mentioned above should not be more than 100%. Therefore, we just introduced this strategy into 10, 30, 50 and 80% coverage rate baselines.

Model fits
The model estimates diagnosed HFMD cases in China for 2015 (n = 1,982,819), 2016 (n = 2,258,846), 2017 (n = 1,948,522) and 2018 (n = 2,398,566), which closely approximates the observed data (Fig. 2). The monthly fitting and observed results are shown in Supplementary  Table 1. The model provides an adequate fit to the observed data as shown in Table 2. All the p values are greater than 0.05, therefore, we cannot reject the null hypothesis at the 5% significant level. The modeled result is a good estimate of the observed value. The model also captures the seasonality of disease incidence with two peaks of HFMD incidences in a given year: March to May and September to October (Fig. 2).

Sensitivity analysis
As we can see from the formula (4), the value of dR 0 dβ is greater than 0 and the values of dR 0 dγ1 , dR 0 dγ2 and dR 0 dγ3 are less than 0. But whether the value of dR 0 dP 1 is greater than 0 remains uncertain. These five parameters all range from 5% below to above of their fitting values. The influence of of β, P 1 and γ 1 on the number of cases are large, while the influence of γ 2 and γ 3 are minimal (Fig. 3).  Table 3 show the impact of bivalent and monovalent vaccines on the incidence of HFMD in one year. Higher vaccine coverage led to reduced incidence but the efficiency also declined as the NNV (number needed to be vaccinated per case reduction) increased. Increased coverage of the monovalent vaccine from 0 to 80% decreased the number of cases by 797, 262(49.1%) in one year. Use of a bivalent vaccine at an 80% coverage level decreased the number of HFMD cases by 828,560(51.0%) in 1 y. Figure 4 shows that when both the EV71 vaccine and the bivalent vaccine reach high coverage, incidence may be fall to zero. As shown in Figs. 4 and 5, a change in bivalent vaccine coverage has a greater impact on I 3 than on I 2 as expected because I 3 reflects the number of cases of HFMD due to both serotypes. Besides, vaccine coverage also has a slight influence on the number of I 2 (Fig. 5). This may indicate that bivalent vaccines are more effective than the sum of their individual effects.
Impact of pulse vaccination strategies Figure 6 and Table 4 show the impact of different strategies of pulse vaccination on the incidence of HFMD. With baseline 10% coverage rate (vaccination rate is 26, 300 person/day), increasing pulse vaccination intensity from 1.0 time to 2.0 times would increase the impact from 130,639 to 167,660 cases prevented, and the NNV would decrease from 73.5 to 71.9. Incidence consistently declined with increasing pulse vaccination intensity across different baseline coverage of the bivalent vaccine. However, the value of NNV showed a slight increase from 92.7 to 93.9 at high baseline coverage of 80%, which is the opposite of other coverage scenarios.

Discussion
There have been many papers on how to control and prevent HFMD from public health or statistical model perspectives. However, there is few work based on mathematical models to simulate dynamics of HFMD and to predict various burdens of HFMD under different scenarios of bivalent vaccine coverage. Our model extends prior work by including vaccination and reproduced the seasonal pattern of HFMD observed among children under 5 years old in China. The findings have important implications for our understanding of the spread of HFMD in China, and consideration of vaccination as a control measure to reduce the spread of HFMD. R 0 is an important indicator of disease outbreak. Its value indicates the number of newly infected persons during one infectious period after one infected person is introduced to an totally susceptible population. If R 0 is larger than 1, the outbreak will take place. If R 0 is under 1, infection will die out. In our study, the formula of R 0 consists of four parts. This is consistent with CC Lai's research [24]. The values of R 0 in China from 2015 to 2018 are all greater than 1. Theoretically it is already over the threshold value for HFMD outbreaks in this period if the virus is introduced and transmitted among the children. There are other studies about the basic reproduction numbers of HFMD, and its numerical range is 1.014-1.742. Wang et al. [25] added the density of pathogen in the contaminated environments as a compartment in their model and they calculated the basic reproduction number R 0 = 1.7424. Li et al. [26] constructed a two-stage structured model to fit the HFMD data from 2009 to 2014 in China and used the same way with us to estimate the basic reproduction number R0. In their study, R 0 fluctuated in the range of 1.06-1.57 in 2019-2014, which is similar to our estimated range.
We also performed a sensitivity analysis to study the effects of different parameters on the incidence of HFMD. Derivative results show that R 0 increases as β increases and γ 1 , γ 2 , γ 3 decrease. Moreover, it's difficult to figure out whether the recovery time courses for people who are hospitalized (1/γ 2 and 1/γ 3 ) are longer than that of people who are not hospitalized (1/γ 1 ). For one, people who are hospitalized may receive effective medical treatment which would shorten the recovery time. For another, people usually go to hospital because of the severity of the condition, which means they may need more time to recover. Therefore, the value of dR 0 dP 1 depends on medical interventions and disease conditions. Considering the biological meanings of these parameters, the trend of R 0 change is consistent with reality. More specifically, the sensitivity analysis curve shows that these five parameters (β, P 1 , γ 1 , γ 2 , γ 3 ) influence the incidence of EV71 and CV-A16 HFMD in different directions, to varying degrees. Considering the biological meanings of β and γ 1 , β, also known as the transmission rate, is composed of two parts: Contact frequency × Infection efficiency. Infection efficiency is determined by the characteristics of different pathogens themselves, and it is difficult for us to intervene. However, we can change the value of contact frequency in a variety of ways. γ 1 indicates the recovery rate of people who are infectious but not hospitalized, which is equal to the reciprocal of recovery time. Most HFMD is a self-limited disease. Therefore, this part of patients who are infectious but not hospitalized accounts for the vast majority of all HFMD patients. Considering that they don't go to hospital, home care and self-purchasing antiviral drugs are key factors influencing the value of γ 1 that we can intervene. Inspiringly, this provides us some new ideas to control HFMD, such as isolation of patients in the early stage of dissemination, more frequent handwashing and ventilation [27] to lower the β value, or better treatment conditions to increase the γ1, γ2 and γ3 value. Apart from constructing the SEIR -VEIR model to fit the reported HFMD data, we introduced vaccination to predict how HFMD incidence changes along with different vaccine coverage rates. It is generally believed that EV71 vaccine can protect against EV71 but not against CVA16 infection, which is also a major cause of HFMD [28]. Therefore, the monovalent vaccine may be less acceptable as a long-term public health strategy to control HFMD. Our findings suggest that bivalent vaccines may be the best strategy to reduce HFMD incidence. Currently, the bivalent vaccine for EV71 and CV-A16 HFMD is under development, the bivalent vaccine will be available in the nearly future. However, the number needed to be vaccinated (NNV) per case reduction suggests that the efficiency of vaccination is low. There are some possible reasons including the low probability of suffering from HFMD among whole population, which is lower than 0.2%. Moreover, the incidence rate in susceptible population (children under 5 years old) is only 2.5% [2]. Secondly, we underestimated the effectiveness of vaccines. Diseases caused by EV71 and CV-A16 are not only HFMD. We did not consider the vaccine protection effect against other diseases.
The pulse vaccination strategy consists of repeated application of vaccine at discrete time with equal interval in a population in contrast to the traditional constant vaccination [20,29]. Pulse vaccination is considered to be a good complement to routine vaccination. Liu [16] integrated the pulse vaccination into a tuberculosis model to vaccinate the susceptible individuals periodically to overcome the shortage of the constant vaccination strategy and the effect of mixed vaccination is obvious, where the proportion of infectious individuals quickly falls down to a very low level. In our study, we  The percentage of the number of HFMD case reduction in I 3 based on the initial number with no intervention ***NNV The number of persons needed to be vaccinated per case reduction introduced pulse vaccination strategy with 6 different intensities, in 4 different baselines in three specific months in view of the apparent seasonal variation of HFMD. There are two peaks of HFMD incidences in one year: May to July and September to October. From the results, we can see the decrease in NNV (number needed to be vaccinated per case reduction) value in 10, 30, and 50% coverage rates along with the pulse intensity increasing, which means as the pulse intensity increases, the efficiency of vaccine also rises up. But in 80% coverage baseline, we find the opposite result, possibly due to the inefficiency of excessive coverage itself.
Although it is somewhat unsatisfying that our analysis could not pinpoint. First, national HFMD epidemic information is calculated monthly. We cannot obtain the daily data which may help us improve the accuracy of the model. Second, serotype replacement indeed exists as the proportions of different HFMD serotypes having been changing in recent years [18], but our model in current study does not account for this transition. We assume that the proportions of different serotypes of HFMD are constant from 2015 to 2018, with reference to baseline levels [19]. Therefore, this may lead to differences from the real HFMD epidemic. Third, EV71 vaccine was launched in China at the end of 2016. But during our investigation, we find people's acceptance of this vaccine is generally low due to its limited effectiveness and its high price. Therefore, some results are just based on our hypothesis but may not be realized in the future considering the cost-effectiveness. Apart from these, we conducted a simple pulse vaccination analysis and just introduced 6 pulse intensities in 3 specific months. However, we have no idea about what intensity of the pulse or pulse vaccination in which months can achieve the maximum of cost-effectiveness.
Overall, as HFMD can be caused by numerous serotypes of enteroviruses, the existing EV71 vaccine has no cross-protection against other subtypes. The epidemic of HFMD is still developing in China. Although EV71 vaccine provides us a more specific way to control the HFMD caused by EV71, which is the leading cause of severe and fatal HFMD, its preventive effect and cost-

2.0×
Case reduction(%*) NNV** Case reduction(%) NNV Case reduction(%) NNV Case reduction(%) NNV Case reduction(%) NNV Case reduction(%) NNV effectiveness on the overall HFMD epidemic is very limited. Other researchers have found that the proportion of EV71-caused HFMD decreased with the popularization of EV71 vaccination, while the proportion of HFMD caused by other serotypes such as CV-A16 and CV-A10 increased [18]. This suggests multivalent vaccine and more effective vaccination strategy different from routine inoculation, such as pulse vaccination, would play an important role in preventing the transmission of all-type-HFMD in the future.