Abstract

Hand, foot, and mouth disease (HFMD), associated with more than 20 disease-causing enteroviruses, is one of the major public health problems in mainland China, and the unique vaccine available is for enterovirus 71 (EV71). In this paper, we propose a new epidemic model to investigate the effect of EV71 vaccination on the spread of HFMD with multiple pathogenic viruses in mainland China. In addition, suitable periodic transmission functions are designed, with a two-year period and taking into consideration the effects of opening and closing of schools. After defining the basic reproduction number , we prove that the disease-free equilibrium is globally asymptotically stable if , and there exists at least one positive periodic solution and the disease is uniformly persistent if . We use the model to simulate the HFMD reported data in mainland China from January 2008 to June 2019. The numerical experiments show that increasing the vaccinated rate can effectively control the spread of HFMD in mainland China, yet the disease does not become extinct. Moreover, if we can control the baseline contact rate of infectious individuals and the recovery rate of symptomatic infectious individuals under certain conditions, which can be achieved by improving protective measures and medical conditions, then the disease will be eliminated.

1. Introduction

Hand, foot, and mouth disease (HFMD) is a common infectious disease among infants and children caused by intestinal viruses of the Picornaviridae family. There exist many types of HFMD virus, such as coxsackievirus A5, A10, A16, A19 types and EV71 type, and so on. Among these viruses, the main viruses causing HFMD are coxsackievirus A16 (CVA16) and enterovirus 71 (EV71) [1]. HFMD mostly occurs in nursery schools or kindergartens, and its high incidence seasons are summer and autumn [1]. The usual incubation period is 2 to 7 days, and the recovery period is 7 to 10 days [2]. Since HFMD was first reported in New Zealand in 1957 [3], it has spread around the world, especially in Asia [47].

HFMD is a mild, self-limiting illness that primarily affects infants and young children. In recent years, HFMD is also a common illness in mainland China. Since 2008, the government has classified it as a class III infectious disease in the National Stationary Notifiable Communicable Diseases (NSNCD) [8], and the monthly cases have been archived by the Chinese Center for Disease Control and Prevention (CCDC) [9] (see Figure 1). The main viruses causing HFMD in mainland China are CVA16 and EV71, and the proportions of HFMD cases caused by CVA16 and EV71 are about 30% and 44%, respectively [1012]. EV71 vaccine, which efficacy against EV71-associated HFMD was 97.4%, has been developed and put into use in mainland China beginning in 2016, so many young children would be protected from EV71 infections [13, 14]. As of May 2019, monovalent EV71 vaccines have not yet been included in the routine paediatric vaccination programme in China, which means that they are needed to be paid out-of-pocket by parents [1517]. Meanwhile, vaccine coverage of these monovalent EV71 vaccines among children aged 6 months to 5 years ranges from to in different provinces in mainland China [1517].

Compartmental differential equation modeling is an important tool for understanding infectious disease spread and control. There are two types of dynamic models for HFMD. One is an autonomous ordinary differential model where all the parameters such as transmission rate, birth population, recovery rate, and so on of whom are assumed as constants. For instance, Tiing and Labadin [2] used a simple model to predict the number of the infected and the duration of an outbreak in Sarawak Malaysia; Roy and Halder [18, 19] proposed a deterministic susceptible-exposed-infectious-recovered () model and a susceptible-exposed-infectious-quarantine-recovered () model of HFMD, and Viriyapong and Wichaino [20] further analyzed dynamic behaviors of ; Yang et al. [21] and Li et al. [22] researched the optimal control of and its dynamic behavior, respectively; Li et al. [23] further constructed a two-stage-structured model to fit the HFMD data in mainland China from 2009 to 2014; Wang et al. [24] further considered indirect transmission coming from the contaminated environments to establish a model ( is the number of the asymptomatic infectious individuals), and Sharma and Samanta [25] considered the combined effect of asymptomatic infectious individuals and quarantine measures on the spread of HFMD. It should be noted that the reported monthly cases of HFMD from CCDC [9] exhibit the periodic pattern. From a practical point of view, clarifying the mechanisms that link seasonal environmental changes to disease dynamics may aid in forecasting the long-term health risks, crucial in developing an effective public health program and in setting objectives, and utilizing limited resources more effectively. For these reasons, the other approach is a nonautonomous differential equation model, which considers periodic transmission rates. Of these, Liu [26] constructed a periodic model to simulate the dynamics of HFMD transmission and showed that quarantine has a positive impact on controlling the spread of HFMD; Zhu et al. [27] used the periodic model to investigate the spread of seasonal HFMD in Wenzhou city, Zhejiang Province, China; Ma et al. [28], considering asymptomatic infectious individuals (), proposed a periodic model to analyse HFMD in Shandong Province, China, which shows that asymptomatic infectious subpopulation plays an important role in the spread of HFMD. However, among the two types of dynamic models, fewer studies considered the influence of EV71 vaccination on the spread of HFMD.

Recently, Samanta [29] and Wang et al. [30] considered vaccination in modeling for HFMD. In these studies, the vaccinated individuals were assumed to be immunized against all HFMD-associated enteroviruses, in which a vaccinated individual would be transferred from the susceptible compartment to the recovered compartment in modeling (see system (1) in [29] or system (1) in [30]). However, in mainland China, only EV71 vaccine existed, and the vaccine only has well immunity for EV71, but not for others [13, 14, 31, 32]. Thus, while there are multiple HFMD-associated enteroviruses, only one vaccine for EV71 is available. A natural problem is how to build a suitable model to reflect the influence of EV71 vaccination on the spread of HFMD with multiple pathogenic viruses in mainland China.

From the practical point of view, periodic transmission functions are needed to consider in modeling. From Figure 1, it should be noted that the monthly reported cases of HFMD in mainland China from CCDC [9] exhibit periodic patterns. The annual cycle patterns or seasonal patterns may be attributed to these facts such as climatic, geographical, and demographic information but not limited to them [28, 30, 33]. Although different seasonality across the whole country, Xing et al.’s epidemiological study of HFMD in China, 2008–2012 [33], implied that geographical differences in seasonal patterns were weakly associated with climate and demographic factors (variance explained 8–23% and 3–19%, respectively). Our analysis of previous studies for periodic HFMD models [4, 26, 28, 30] finds that all the chosen periodic transmission functions have a similar form, and its period is 1 year, which effectively can reflect the annual seasonal changes for HFMD. However, annual cycle patterns aside, the epidemics of HFMD in mainland China have displayed the following two phenomena: (1) the data from CCDC [8] (see Figure 1) show that HFMD has a two-year epidemic cycle because the epidemic period of partial HFMD-associated enteroviruses in some areas is two years [34, 35]; (2) Figure 1 shows that the number of cases suddenly rises in September and October each year because of the opening and closing of schools in mainland China [36]. The previous models based on the disease’s transmission functions [4, 26, 28, 30] cannot simulate the above two phenomena. Hence, to better reveal the mechanism of the epidemic, it is very important to find a more suitable periodic transmission function for modeling.

Given the issues discussed above, we propose a new mathematical model to investigate the effect of EV71 vaccination on the spread of HFMD in mainland China. To analyze the effect of EV71 vaccination, we derive our modeling ideas from the previous modeling for influenza [37]. Then, in modeling, we divide the total individuals into two groups, the unvaccinated and the vaccinated. In addition, novel two-year periodic transmission functions, with considering the effects of school opening and closing, are also considered in our modeling. Then, we evaluate the basic reproduction number to analyze the dynamical behaviors of the model and use the proposed model to fit the monthly reported data of HFMD in mainland China from January 2008 to June 2019. Furthermore, through sensitivity analysis of the basic reproduction number in terms of key parameters, we explore some effective prevention and control measures for HFMD in mainland China.

The paper is organized as follows: in Section 2, we introduce a new periodic epidemiological model of HFMD. In Section 3, the dynamic behaviors of this model are analyzed theoretically. In Section 4, simulations of the model, sensitivity analysis of the basic reproduction number, and some prevention and control measures are performed. In Section 5, we discuss and summarize our conclusions.

2. Model Formulation

Recently, Samanta [29], Wang et al. [30] considered the effect of the vaccination in modeling for HFMD, in which it assumed that vaccine can resist all the HFMD-associated enteroviruses. However, there only exists EV71 vaccine to be used in mainland China, and the vaccine only has well immunity for EV71 but not for others [13, 14]. It notes that the modeling techniques in [29, 30] cannot be used to really reflect the spread of the disease with a single effective vaccine and multiple viruses. In order to model the influence of EV71 vaccination on the spread of HFMD with multiple pathogenic viruses in mainland China, we divide the total individuals into two different groups, the unvaccinated and the vaccinated, in modeling. Next, based on the above analyses, we denote the total number of unvaccinated and vaccinated individuals by and , respectively, and classify each of them into four subclasses: susceptible, exposed, infectious with symptoms (symptomatic infectious), and infectious but not yet symptomatic (asymptomatic infectious), with the number of unvaccinated individuals denoted by and and the number of vaccinated individuals denoted by and , respectively. Denote the total population size .

In addition, many epidemiological models for HFMD ([2628, 30]), considering annual seasonal changes, were simulated by using the sinusoidal function with one year period (take a month as the unit time). As discussion in Section 1, seasonal changes and school holidays will change the contact rate of children, and the virus contagion changes periodically, some of which are one year for a cycle and some of which are two years for a cycle. Thus, taking a month as the unit time, we define periodic transmission functions with the period of 2 years as and for symptomatic infectious individuals and asymptomatic infectious individuals, respectively. The studies in [31] showed that, in mainland China, the EV71 vaccine was targeted at infants or children aged 6 to 71 months. They were also a high-risk group for HFMD [31, 32]. Next, let be the monthly number of children entering the age of 6–71 months, and d be progression rate leaving the children group aged 6–71 months. We also assume that p is the vaccinated rate for infants or children aged 6–71 months. Moreover, the results in [15, 17] showed that vaccine coverage of these monovalent EV71 vaccines among children aged 6 months to 5 years ranges from to in different provinces in mainland China in recent years. This implies that the current vaccination rate of the monovalent EV71 vaccine is very low in mainland China. In the case of the low vaccine coverage, we let be the proportion of HFMD cases in mainland China, without vaccination, caused by HFMD-associated enteroviruses except EV71. Furthermore, we assume that a symptomatic or an asymptomatic infectious individual with vaccination has the same recovery rate as a symptomatic or asymptomatic infectious individual without vaccination. Then, the transmission dynamics associated with these subpopulations are shown in Figure 2.

The model is described by the following ordinary differential equations:where p is nonnegative and other parameters are positive, and its biological meanings are listed in Table 1.

3. Stability Analysis and Persistence

In this section, we will investigate the dynamic behaviors of system (1). First, we introduce some notations which will be used throughout this paper. denotes the n-dimensional Euclidean space, and . denotes a real matrix. The superscript T denotes matrix or vector transposition. Letbe the initial condition of system (1), andbe the solution of system (1) at time t for .

From the first equation of system (1), if , then . According to function continuity and monotonicity, when , we have . Moreover, we have similar results for other state variables of system (1). Therefore, any solution of system (1) with nonnegative initial conditions is nonnegative. Next, the following lemma shows that the solutions are uniformly ultimately bounded.

Lemma 1. The solutions of system (1) eventually enterMoreover, is a positively invariant set.

Proof. It is obvious that the unvaccinated population size and the vaccinated population size . From system (1), we haveAccording to the comparison theorem, it follows from (5) and (6) thatUsing (7) and (8), we obtain that and . Hence, the solutions of system (1) are uniformly ultimately bounded. Moreover, if the initial condition , then the solution , i.e., is a positively invariant set for system (1). This completes the proof.
It is obvious that system (1) always has a disease-free equilibrium:Next, we first introduce a very important threshold basic reproduction number . The basic reproduction number is the number of secondary cases in which one case would produce in a completely susceptible population [40]. We calculate the basic reproduction number for system (1) by using the general calculation procedure given by Wang and Zhao [41]. It is easy to verify that system (1) satisfies the conditions (A1)–(A7) given in [41], and we haveThus, we obtain thatLet be the monodromy matrix of the linear ω-periodic system , and be the spectral radius of , where ω is the period. Let be the matrix solution of the following initial value problem:where I is the identity matrix.
Suppose that is the initial distribution of infectious individuals, then is the rate of new infections produced by the infectious individuals introduced at time s, and represents the distribution of those infectious individuals who are newly infected at time s and still remain in the infected compartment at time t for . Then, one has thatis the distribution of accumulative new infections at time t produced by all those infected individuals introduced at time previous to t. Let be the ordered Banach space of all ω periodic functions from to with maximum norm and positive cone . Then, we can introduce a linear operator byFollowing the results obtained in [41], the basic reproduction number for model (1) is defined as the spectral radius of operator L, i.e., .
In order to analyze the dynamic behaviors of system (1), we have the following result.

Lemma 2 (Theorem 2.2 in [41]). The following statements are valid:(i) if and only if (ii) if and only if (iii) if and only if Therefore, the disease-free equilibrium of system (1) is locally asymptotically stable if and unstable if . Next, we will analyze the global dynamics of the disease-free equilibrium .

Theorem 1. If , then the disease-free equilibrium of system (1) is globally asymptotically stable.

Proof. According to Lemma 2, one has that is equivalent to . Hence, we only need to prove that when , is globally asymptotically stable.
Choose small enough such that , whereFrom Lemma 1, for , there exists such that for , and . Thus, for , it follows from system (1) thatConsider the following auxiliary system:where . Following Lemma 2.1 in [42], there exists a positive ω-periodic function such that is a solution of system (17), where . Due to the continuity of functions and following Lemma 2, when and , then . Hence, we have as , which implies that the zero solution of system (17) is globally asymptotically stable. Based on the comparison principle, it follows from system (17) thatUsing (18), it follows from the fifth and last equations in system (1) thatUsing (18) and (19), it follows from the first and sixth equations in system (1) thatThus, is globally asymptotically stable when . This completes the proof.
Next, we will analyze the persistence of system (1). DefineAssume that is the solution of system (1) with the initial value . According to the fundamental existence-uniqueness theorem in [43], is unique. Denote that is the Poincaré map with respect to system (1), i.e.,where ω is the period. It is obvious thatIt follows from Lemma 1 that is positively invariant and P is point dissipative for system (1). In order to analyze the persistence for system (1), we give the following lemma.

Lemma 3. If , then there exists such that for any with , we havewhere D is a distance function in .

Proof. According to Lemma 1, if , then . Choose small enough such that . Now, we prove thatIf it is false, thenfor some . Assume that there exists such that . According to the continuity of the solutions associated with the initial conditions, when , it follows thatFor , let , where and , which is the greatest integer less than or equal to . Then, for , we havewhich implies that . Then, for , one hasConsider the following auxiliary linear system:As the similar proof process in Theorem 1, there exists a positive ω-periodic function such that is a solution of system (30), where . Since , it follows that if , then as . From the comparison principle, we haveasfor , which is a contradiction with Lemma 1. Therefore, . This completes the proof.

Theorem 2. If , then system (1) is uniformly persistent and has at least one positive periodic solution.

Proof. First, we prove that system (1) is uniformly persistent, i.e., there exists such that any solution of system (1) with the initial value satisfiesFor , solving system (1), we obtain thatFrom (34), it is obvious that for , the solution for system (1) as the initial value . This means that is positively invariant. In addition, is relatively closed in , and following Lemma 1, the discrete-time system admits a global attractor in . Then, we setNow, we prove thatIt is easy to know that , then we only need to prove . That is, for anywe haveIf not, then there exists such thatIf we replace the initial time with , then it follows from (34) thatfor . This means that , for , which contradicts the definition of . Therefore, it follows that . Moreover, there only exists one fixed point of P in .
According to Lemma 3, is an isolated invariant set in . From the acyclicity theorem on uniform persistence for maps (Theorem 1.3.1 in [44]), it follows that P is uniformly persistent with respect to . Moreover, according to Theorem 1.3.1 in [44], the solutions of system (1) are uniformly persistent with respect to , i.e., there exists such that any solution of system (1) with the initial value satisfiesNext, we prove that system (1) has a positive periodic solution. Theorem 1.3.6 in [44] implies that P has a fixed pointMoreover, we can also obtain that for system (1). Suppose not, assume . By the periodicity of P, we have . From the first equation of (34), one haswhich contradicts with . Thus, it follows that . Then, from (34), it is easy to obtain that , for . Due to the definition of semiflow P, we obtain thatis a positive ω-periodic solution of system (1). This completes the proof.

4. Simulations and Sensitivity Analysis

In this section, by using system (1), we simulate the reported data of HFMD in mainland China from January 2008 to June 2019, predict the trend of the disease, and seek for some prevention and control strategies. We obtain the monthly number of newly reported HFMD cases in mainland China from the website of the Chinese Center for Disease Control and Prevention (CCDC) [9] and the National Health and Family Planning Commission of the People’s Republic of China (NHFPC) [8]. We only consider the population of young children aged 6–71 months, who are a high-risk group for HFMD [31, 32] and get demographic information from the National Bureau of Statistics of China (NBSC) [38]. From NBSC [38], we obtain the monthly number of children entering the age of 6–71 months and the average monthly rate of leaving the children group aged 6–71 months . From CCDC [9] and NHFPC [8], we get the disease-related death rate . The proportion of infectious individuals without vaccination caused by HFMD viruses except EV71 is [39]. Other parameters can be obtained from the literature or assumed on the basis of common sense (see Table 1). Unfortunately, there is no public data on the vaccinated rate p and the proportion of symptomatic infectious individuals ρ. Therefore, we estimate the values of unknown parameters containing and . By using the BP neural network algorithm [45] (see Appendix), we estimate parameters of system (1) by calculating the following error function:where is the reported data of HFMD at time t. is the numerically computed solution of of system (1) at time t. n is the total number of sample points collected at all times.

Next, we set the initial conditions for system (1). We take January 2008 as the initial time. It notes that the EV71 vaccine has been vaccinated in mainland China starting in 2016, then vaccinated rate p of system (1) is zero before 2016 [13, 14]. Hence, since the total number of vaccinated individuals before January 2016, then it follows that . The initial number of the susceptible human population aged 6–71 months of 2008 is calculated as from NBSC [8], and the number of the initial unvaccinated symptomatic infectious humans is from CCDC [9]. However, the initial number of the unvaccinated exposed humans , the unvaccinated asymptomatic infectious humans , and the unvaccinated recovered humans cannot be obtained. We derive reversely by the parameter ; get by using the proportion parameter ρ to compare with , and is estimated roughly. The numerical simulation of the model on the number of HFMD cases in mainland China from January 2008 to June 2019 is shown in Figure 3, and its fitting error (45) is calculated as . In addition, in order to test how well system (1) actually reflects the data, by using the method in [22, 46], we consider the following hypotheses.

Null hypothesis, : the estimated parameters of system (1) are equal to actual values.

Alternative hypothesis, : the estimated parameters of system (1) are not equal to actual values.

Since the EV71 vaccine was only started in 2016, we consider hypothesis test for the two time periods from January 2008 to December 2015 and January 2016 to June 2019, and the chi-square values and degrees of freedom are shown in Table 2. Then, from Table 2, we obtain that it cannot reject the null hypothesis at the 5% significant level by Pearson’s criterion of chi-square test [47]. Therefore, it indicates that, with these parameter values, there is a good fit between the simulation of model (1) and the HFMD cases in mainland China. Furthermore, based on our parameters listed in Table 1, we calculate the basic reproduction number as after 2015, for this case where the vaccination is considered, which means that system (1) admits at least a positive periodic solution (see Theorem 2).

It notes the initial conditions except adopted in model fitting are mostly assumed and back-extrapolated by parameters. Thus, it is necessary to analyze the influence of the initial values on the number of infected cases , which is illustrated in Figure 4. Figure 4 shows that the initial value has a greater impact on while other initial values have little or no impact on . This implies that the selection of the initial values is very important for our simulation. According to [6, 7, 10, 11, 32], due to most of the infected individuals of HFMD are 6–71 months old, then we carefully select that the susceptible individuals are 6–71 months old.

In order to investigate the effect of the EV71 vaccination strategy on HFMD in mainland China, we compare the trends of HFMD between with vaccination and without vaccination by simulating model (1), which is presented in Figure 5. In Figure 5, the solution curve without considering vaccination is on top of the one with considering vaccination after 2015, which means that the EV71 vaccination strategy is effective in reducing the spread of HFMD in mainland China in the last two years. In addition, we calculate the basic reproduction number of system (1) without considering vaccination () as 1.24 and with considering vaccination () as 1.09, so system (1) in each case admits a positive periodic solution according to Theorem 2 (see the solution curves from 2030 to 2035 in Figure 5). In addition, Figure 6 shows that under the current vaccination measures, the total number of patients in 2016, 2017, and 2018 has reduced about , , and , respectively.

To compare the performance of the periodic transmission functions in this paper with some previous research studies, by using system (1), we consider the following two cases:

Periodic functions (I): choose the transmission functions in this paper such as and , where the values of are given in Table 1.

Periodic functions (II): choose the transmission functions in the previous research studies [2628, 30] such as and . Based on our parameter values in Table 1, we estimate that the optimal values of are , respectively.

The studies in [48] showed that the trend of monthly new cases of infectious diseases is usually closely related to the transmission functions. It notes that the periods of and are 12 months, and these functions have only one peak in each year. Meanwhile, the periods of and are 24 months, and these functions have two peaks in each year. As explained in Introduction, the transmission functions are closer to the trend of actual data change in mainland China than . Moreover, in Figure 7, the simulation of of system (1) is based on the periodic functions (I) (i.e., ) but the periodic functions (II) (i.e., ) can well fit the two phenomena, where HFMD has a two-year epidemic cycle and the number of the cases per year has a slight upward trend in September and October in mainland China (see the reasons in Introduction).

Moreover, the period of the solution for system (1) with periodic functions (II) is 12 months, which fails to reflect the 24 months periodic variation of the disease, while these biological phenomena are well presented by the solution for system (1) with periodic functions (I). Furthermore, we calculate the fitting error for system (1) with periodic functions (I) and with periodic functions (II) which is and , respectively. On the whole, our model, especially based on our selection of the transmission rate functions, can well reveal the objective laws of the spread of HFMD in mainland China.

According to Theorems 1 and 2, we have known that the basic reproduction number , which determines the outbreak of the disease or extinction, is an important quantity in characterizing the spread of disease. Therefore, to find some effective prevention and control measures for HFMD in mainland China, we perform sensitivity analysis to determine the influence of some parameters on . Due to some parameters of system (1) for HFMD are difficult or not to control such as and so on, biologically, then we only consider the influence of the parameters on , which are depicted in Figure 8, and its related prevention and control measures are shown in Figure 9. Under the assumption of low vaccine coverage, choosing , Figure 8(a) shows that decreases with the increase of p. This implies that the larger vaccinated rate p () is, the less number of symptomatic infectious individuals is (see Figure 9(a)). Moreover, Figure 8(a) also shows that is always larger than 1 as , which means that without other prevention and control, the disease is not extinct. Figures 8(b) and 8(c) show that the less the baseline contact rate of or is, the less is. Moreover, when the baseline contact rate or , can be less than 1, which means that the disease is extinct. However, based on our estimating parameter values, or is far greater than or , which means that HFMD becomes an endemic disease in mainland China. If the relevant departments can effectively control the contact rate between the healthy children and the infectious children and improve hygienic precautions and environmental cleaning, such as washing hands before meals and after using the toilet, and making air fresh indoors and so on, the disease can be effectively control or extinct (see Figures 9(b) and 9(c)). Figure 8(d) shows that decreases with the increasing of the recovery rate of symptomatic infectious individuals , i.e., if average recovery period of symptomatic infectious individuals is less, then is less, and the disease is extinct if . Figure 9(d) shows that the number of symptomatic infectious cases decreases as the recovery rate increases, especially the disease is extinct when . Therefore, we suggest that the patients should seek medical advice in time to achieve diagnosis and treatment early. At the same time, because the disease does not have targeted drugs, some targeted drugs need to be developed actively. Finally, in order to further discuss the credibility of our results, we will evaluate and discuss the consistency of the scope of with other studies in China. Although the ideas of dynamic modeling for different studies based on the HFMD cases in China are not the same, the range of of models may have similar characteristics. For instance, the study in [30], which assumed that the vaccine was against all HFMD-associated enteroviruses and analyzed the HFMD cases in mainland China from 2010 to 2014, showed that as the vaccination rate, the baseline contact rate of symptomatic infectious individuals and the baseline contact rate of asymptomatic infectious individuals change, the ranges of are , , and , respectively. The study in [28], based on the HFMD cases in Shandong Province, China, which mainly discussed the impacts of quarantine measures and the fraction of asymptomatic infectious individuals on the spread of HFMD, indicated that, as the quarantine rate and the fraction of asymptomatic infectious individuals change, the ranges of are and . In a word, most studies for HFMD in China (see paper [27, 28, 30, 49]) showed that the range of is around 1 and not far away from 1, which is similar to the change of shown in Figure 8. Therefore, these discussions can also indirectly improve the credibility of our results.

5. Discussion

In this paper, we have proposed and analyzed a HFMD model with periodic transmission rates to take into account seasonal outbreak of HFMD infection and to consider EV71 vaccination. The studies in [29, 30] have assumed that the vaccine is effective against all enteroviruses, which may not be in line with the current situation of HFMD transmission in mainland China. In fact, one vaccine is available in mainland China and it works only for EV71 but not for others. To investigate the effects of EV71 vaccination on HFMD in the presence of multiple pathogenic viruses, our proposed model divided the total individuals into the unvaccinated and the vaccinated, which were described in different compartments (see Section 2). Moreover, we have chosen novel periodic transmission functions, unlike the previous studies in [4, 26, 28, 30], to reflect the two-year periodic nature and the infections’ change with the opening and closing of schools in mainland China (see Sections 1 and 2).

We have analyzed the global dynamics of system (1) by evaluating the basic reproduction number : if , then the disease-free equilibrium of system (1) is globally asymptotically stable (see Theorem 1), which means that the disease is extinct; if , then system (1) uniformly persists and has at least one positive periodic solution (see Theorem 2), which means that the disease persists or has not been eliminated.

By fitting our model to the reported data on symptomatic cases of HFMD in mainland China, we estimated the basic reproduction number from 2008 to 2015 and from 2016 to June 2019. The major factors for this difference are that EV71 vaccination has been used since 2016, the reason of which can also be seen in the analyses for the influence of the vaccinated rate p on in Figure 8(a). Moreover, the periodic transmission functions chosen in our paper perform efficiently for system (1) to fit the reported data in mainland China, given the analyzes of Figure 7. Furthermore, to find effective prevention and control measures for HFMD in mainland China, we did sensitivity analysis to determine the influence of some artificial control parameters on . It showed that under the assumption of low vaccine coverage, increasing the vaccinated rate with respect to EV71 vaccination can reduce the number of infections, but the disease does not go extinct. If we can effectively control the contact rate between the healthy children and the infectious children (see the analyses of Figures 8(b) and 8(c) and Figures 9(b) and 9(c)) and improve the medical conditions (see the analyses of Figures 8(d) and 9(d)), the disease may be eliminated. In a word, the proposed prevention and control measures may be useful for public health governance.

It notes that the observed aggregated monthly HFMD cases in Figure 1 are a mixture of different seasonality in mainland China. Also, Xing et al.’s study of HFMD in mainland China from 2008 to 2012 [33] showed that seasonal diversities (although this was not the only factor) made the spread of HFMD present different characteristics in different regions, such as the annual amplitude of HFMD epidemics increased with increasing latitude and semiannual periodicity was the strongest in the south, and so on. It is difficult for us to use the data in Figure 1 to discuss the differences in the spread of HFMD in different regions of the country. Thus, we used the national average transmission rate functions in our modeling, where are periodic functions. In future studies, in order to avoid the influence of seasonal differences in different regions on the research results, we will consider a specific region such as a province in China.

Appendix

The BP Neural Network Algorithm for Estimating Parameters

Based on the BP neural network, we design an algorithm to estimate unknown parameters of system (1). Denote the unknown parameter vector for system (1) (see Table 1). Let and (the number 138 is the total months from January 2008 to June 2019). Then, the estimation algorithm is described as follows (the detailed process of the BP neural network can be seen in [45], and the topology of the BP neural network is shown in Figure 10):

Step 1: generate sample data: generate n group sample data , by Latin hypercube sampling. Then, taking into system (1), we obtain n group numerically computed data , (see Figure 11).Step 2: use the BP neural network to train data: set and , and each element of and is normalized into interval . Then, taking as the input vector and as the output vector, the BP neural network is used to train dataset, and let be the training completed network, which the training error is less than the setting value (see Figure 12).Step 3: estimate parameters: first, taking reported data of HFMD as the input vector into the training completed network , we gain an output forecast vector . Then, taking into system (1), it gets the simulation of . Finally, we set a small positive content to estimate the error function (45). If , the estimated parameter vector is reasonable; otherwise, it goes to step 1 and reset the settings of the BP neural network.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was partially supported by the National Natural Science Foundation of China (nos 11571170 and 11971013) and the Joint Special Funds for Basic Research in Local Undergraduate Universities (Part) of Yunnan Province of China (nos 2018FH001-113).