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

Abstract

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 R₀ 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.

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 laboratory-confirmed 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₂); infectious hospitalized with EV71-HFMD or CVA16-HFMD (I₃) 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₂ and I₃ compartments in 2017 and 2018, where ${\mathrm{I}}_{\mathrm{h}}^{-\mathrm{e}}$ indicates people who are infectious with pathogens other than EV71 and ${\mathrm{I}}_{\mathrm{h}}^{\mathrm{e}}$ 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₂ + I₃ + R. Fig. 1 and Table 1 describe the model parameters.

Fig. 1

An extended SEIR model of Hand, foot and mouth diseases including 6 or 8 compartments. Figure 1A indicates the dynamic of HFMD transmission under no vaccination; Fig. 1B indicates the dynamic of HFMD transmission under vaccination

Table 1 Model parameters, description, source and the basic reproduction number in China, 2015–2018

The dynamics of HFMD is described in the following set of coupled ordinary differential equations to represent the dynamics of HFMD without vaccination.

$$ \left\{\begin{array}{c}\ \frac{dS}{dt}=\pi N+\gamma R-\frac{\beta S\sum I}{N}-\mu S\ \\ {}\ \frac{d{E}_s}{dt}=\frac{\beta S\sum I}{N}-\alpha \left({P}_1+{P}_2+{P}_3\right){E}_s-\mu {E}_s\ \\ {}\ \frac{d{I}_N}{dt}=\alpha {P}_1{E}_s-{\gamma}_1{I}_N-\mu {I}_N\ \\ {}\ \frac{d{I}_2}{dt}=\alpha {P}_2{E}_s-{\gamma}_2{I}_2-\mu {I}_2\\ {}\ \frac{d{I}_3}{dt}=\alpha {P}_3{E}_S-{\gamma}_3{I}_3-\mu {I}_3\ \\ {}\ \frac{dR}{dt}={\gamma}_1{I}_N+{\gamma}_2{I}_2+{\gamma}_3{I}_3-\gamma R-\mu R\end{array}\right. $$

(1)

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).

$$ \left\{\begin{array}{c}\frac{dS}{dt}=\pi N+\gamma R-\frac{\beta S\sum I}{N}-v-\mu S\\ {}\frac{dV}{dt}=v-\frac{\beta V\sum I}{N}-\mu V\\ {}\frac{d{E}_s}{dt}=\frac{\beta S\sum I}{N}-\alpha \left({P}_1+{P}_2+{P}_3\right){E}_s-\mu {E}_s\\ {}\frac{d{E}_v}{dt}=\frac{\beta V\sum I}{N}-\alpha \times \left[{P}_1+{P}_2+\left(1- VE\right)\times {P}_3\right]{E}_v-\mu {E}_v\ \\ {}\frac{d{I}_N}{dt}=\alpha {P}_1\left({E}_S+{E}_v\right)-{\gamma}_1{I}_N-\mu {I}_N\\ {}\frac{d{I}_2}{dt}=\alpha {P}_2\left({E}_S+{E}_v\right)-{\gamma}_2{I}_2-\mu {I}_2\ \\ {}\frac{d{I}_3}{dt}=\alpha {P}_3{E}_S+\left(1- VE\right)\alpha {P}_3{E}_v-{\gamma}_3{I}_3-\mu {I}_3\\ {}\frac{dR}{dt}={\gamma}_1{I}_N+{\gamma}_2{I}_2+{\gamma}_3{I}_3-\gamma R-\mu R\end{array}\right. $$

(2)

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₂(0) + I₃(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 age-groups. 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=\sum {\left({\mathit{\log}}_2\left(1+ data\ per\ month\right)-{\mathit{\log}}_2\left(1+{I}_2\ and\ {I}_3\ on\ day1+ on\ day2+\dots + on\ day28, or29,30,31\right)\right)}^2 $$

(3)

Where I₂ indicates the number of no EV71 or CVA16 HFMD $\Big({\mathrm{I}}_{\mathrm{h}}^{-\mathrm{ec}}$) / no EV71 HFMD $\left({\mathrm{I}}_{\mathrm{h}}^{-\mathrm{e}}\right)$, I₃ indicates the number of EV71 or CVA16 HFMD $\left({\mathrm{I}}_{\mathrm{h}}^{\mathrm{ec}}\right)$ / only EV71- HFMD $\left({\mathrm{I}}_{\mathrm{h}}^{\mathrm{e}}\right)$.

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₀: The modeled results (I₂ + I₃) are equal to the observed number of HFMD cases shown in Table S1. Alternative hypothesis H₁: The modeled results are not equal to the observed values shown in Table S1.

Derivation of basic reproductive number (R₀)

The basic reproductive number (R₀) 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₀) 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₀ is related to transmission rate (β), progression rate (α), the proportion of I_N, I₂, I₃ (P₁, P₂, P₃), the recovery rate (γ₁, γ₂, γ₃), the mortality rate (μ) and the vaccine efficacy (VE).

Sensitivity analysis

Considering the formula of R₀(Formula (2) in the Additional file 2, the expression of R₀ in 2015 and 2016), the basic reproductive number is affected by β, α, P₁, γ₁, γ₂ and γ₃ while α is fixed to 1. Therefore, we take a derivation with respect to β, P₁, γ₁, γ₂ and γ₃, respectively. Derivation results are as follows:

$$ \left\{\begin{array}{c}\frac{d{R}_0}{d\beta}=\frac{1}{\alpha +\mu}\times \left(\frac{\alpha {P}_1}{\gamma_1+\mu }+\frac{0.3\times \alpha \left(1-{P}_1\right)}{\gamma_2+\mu }+\frac{0.7\times \alpha \left(1-{P}_1\right)}{\gamma_3+\mu}\right)\\ {}\frac{d{R}_0}{d{P}_1}=\frac{\beta \alpha}{\alpha +\mu}\times \left(\frac{1}{\gamma_1+\mu }-\frac{0.3}{\gamma_2+\mu }-\frac{0.7}{\gamma_3+\mu}\right)\ \\ {}\frac{d{R}_0}{d{\gamma}_1}=-\frac{\beta \alpha {P}_1}{\left(\alpha +\mu \right)\times {\left({\gamma}_1+\mu \right)}^2}\\ {}\frac{d{R}_0}{d{\gamma}_2}=-\frac{0.3\times \beta \alpha \left(1-{P}_1\right)}{\left(\alpha +\mu \right)\times {\left({\gamma}_2+\mu \right)}^2}\ \\ {}\frac{d{R}_0}{d{\gamma}_3}=-\frac{0.7\times \beta \alpha \left(1-{P}_1\right)}{\left(\alpha +\mu \right)\times {\left({\gamma}_3+\mu \right)}^2}\ \end{array}\right. $$

(4)

We further conducted a sensitivity analysis to measure the impact of these five parameters (β, P₁, γ₁, γ₂ and γ₃) on the number of EV71 or CV-A16 HFMD (I₃) 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₂ and I₃ under different vaccine coverage scenarios. To compare HFMD incidence under the monovalent EV-71 vaccine and bivalent vaccine coverage scenarios, we compared scenarios via a 5 × 6 vaccination matrix as follows: a) No EV-71 vaccination throughout the year; b) 10% EV-71 vaccination coverage (vaccination rate: 26,300 person/day); c) 30% EV-71 vaccination coverage (vaccination rate: 78,900 person/day); d) 50% EV-71 vaccination coverage (vaccination rate: 131,500 person/day); e) 80% EV-71 vaccination coverage (vaccination rate: 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.

Results

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).

Fig. 2

Fitting result between the simulated number and the monthly reported cases of hand, foot and mouth disease in China, 2015–2018. The observed number of HFMD cases comes from official infectious disease report from China CDC, the fitting number of HFMD cases is the sum of I₂ and I₃, where I₂ indicates the number of no EV71 or CVA16 HFMD $\Big({\mathrm{I}}_{\mathrm{h}}^{-\mathrm{ec}}$) / no EV71 HFMD $\left({\mathrm{I}}_{\mathrm{h}}^{-\mathrm{e}}\right)$, I₃ indicates the number of EV71 or CVA16 HFMD $\left({\mathrm{I}}_{\mathrm{h}}^{\mathrm{ec}}\right)$ / only EV71- HFMD $\left({\mathrm{I}}_{\mathrm{h}}^{\mathrm{e}}\right)$

Table 2 Chi-square test results of fitting model, 2015–2018

Sensitivity analysis

As we can see from the formula (4), the value of $\frac{d{R}_0}{d\beta}$ is greater than 0 and the values of $\frac{d{R}_0}{d\gamma 1}$, $\frac{d{R}_0}{d\gamma 2}$ and $\frac{d{R}_0}{d\gamma 3}$ are less than 0. But whether the value of $\frac{d{R}_0}{d{P}_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₁ and γ₁ on the number of cases are large, while the influence of γ₂ and γ₃ are minimal (Fig. 3).

Fig. 3

Estimated impacts of 5 parameters on the number of EV71 or CV-A16 HFMD cases (I₃), with ranges based on sensitivity analysis. A, B, C, D and E: The five parameters (β, P₁, γ₁, γ₂, γ₃) range from 5% below to above of their fitting values, respectively

Impact of different vaccination type and coverage

Figures 4 and 5 and 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₃ than on I₂ as expected because I₃ reflects the number of cases of HFMD due to both serotypes. Besides, vaccine coverage also has a slight influence on the number of I₂ (Fig. 5). This may indicate that bivalent vaccines are more effective than the sum of their individual effects.

Fig. 4

Number of HFMD cases in I₃ changes under different EV71 and bivalent vaccine coverage rates. A-E: 0, 10, 30, 50, 80% EV71 vaccine coverage combined with 0, 10, 30, 50, 80, 100% Bivalent vaccine coverage, respectively. I₃ indicates the number of EV71 or CVA16 HFMD $\left({\mathrm{I}}_{\mathrm{h}}^{\mathrm{ec}}\right)$ / only EV71- HFMD $\left({\mathrm{I}}_{\mathrm{h}}^{\mathrm{e}}\right)$, based on different types of vaccines

Fig. 5

Number of HFMD cases in I₂ changes under different EV71 and bivalent vaccine coverage rates. A-E: 0, 10, 30, 50, 80% EV71 vaccine coverage combined with 0, 10, 30, 50, 80, 100% Bivalent vaccine coverage, respectively. I₂ indicates the number of no EV71 or CVA16 HFMD $\Big({\mathrm{I}}_{\mathrm{h}}^{-\mathrm{ec}}$) / no EV71 HFMD $\left({\mathrm{I}}_{\mathrm{h}}^{-\mathrm{e}}\right)$, based on different types of vaccine

Table 3 Total number of HFMD case reduction in I₃* under different coverage scenarios over 1 year

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.

Fig. 6

Changes of EV71 or CV-A16 HFMD (I₃) cases under different pulse vaccination scenarios. A, B, C, D: 6 pulse intensities (1.0×, 1.2×, 1.4×, 1.6×, 1.8×, 2.0×) based on 4 different vaccination rates (26,300 person/day, 78,900 person/day, 131,500 person/day and 210,400 person/day), respectively

Table 4 Numbers of EV71 or CV-A16 HFMD case reduction under different pulse vaccination scenarios over 1 year

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₀ 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₀ is larger than 1, the outbreak will take place. If R₀ is under 1, infection will die out. In our study, the formula of R₀ consists of four parts. This is consistent with CC Lai’s research [24]. The values of R₀ 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₀ = 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₀ 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₀ increases as β increases and γ₁, γ₂, γ₃ decrease. Moreover, it’s difficult to figure out whether the recovery time courses for people who are hospitalized (1/γ₂ and 1/γ₃) are longer than that of people who are not hospitalized (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 $\frac{d{R}_0}{d{P}_1}$ depends on medical interventions and disease conditions. Considering the biological meanings of these parameters, the trend of R₀ change is consistent with reality. More specifically, the sensitivity analysis curve shows that these five parameters (β, P₁, γ₁, γ₂, γ₃) influence the incidence of EV71 and CV-A16 HFMD in different directions, to varying degrees. Considering the biological meanings of β and γ₁, β, 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. γ₁ 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 γ₁ 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 hand-washing 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 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-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.