Abstract

In this paper, we investigate the dynamics and optimal control strategies of a modified hand, foot, and mouth disease (HFMD) model incorporating the EV-A71 vaccination in Wenzhou, China, analytically and numerically. We define the basic reproduction number and show that it can be used to determine whether HFMD becomes extinct or not. Based on the monthly reported HFMD cases in Wenzhou for 76 months, we estimate the parameters in the dynamic model by using the method of minimum chi-square fitting, conduct the sensitivity analysis to investigate the influence of each uncertain parameter on with the methods of Latin hypercube sampling and partial rank correlation coefficient, and find that the EV-A71 vaccination does not lead to the extinction of HFMD, but slightly reduces the incidence of HFMD. In order to control the spread of HFMD in Wenzhou, we need to increase the rate of EV-A71 vaccination, decrease the contact rates, and shorten the course of disease.

1. Introduction

Hand, foot, and mouth disease is a common infectious disorder caused by various enteroviruses, and enterovirus 71 (EV71) and Coxsackievirus A16 (CV-A16) are the most commonly reported [13]. HFMD predominantly affects children younger than 5 years old, and most patients exhibit a self-limiting illness that typically includes fever, skin eruptions on the hands and feet, and vesicles in the mouth [1, 2, 4]. Also, it is a highly contagious disease with a latency period of 2–7 d. More often than not, the patients will recover over 7 to 10 d [5, 6].

There are no specific treatment drugs and vaccine available for HFMD, but three inactivated monovalent EV-A71 vaccinations have been licensed in mainland China, all of which have demonstrated high efficacy (90.0%–97.4%) against EV-A71-associated HFMD in infants and young children. However, these vaccinations do not offer protection against HFMD caused by the CV-A16 serotype or others [2, 7, 8].

In recent years, many dynamical HFMD models have been revealed as a powerful tool to analyze the spread and control of HFMD qualitatively and quantitatively [915]. Ma et al. [15] and Li et al. [16, 17] pointed out that numerous subclinical cases (adults) who carry HFMD virus but have no symptoms play an important role, leading to the recurrent outbreaks of HFMD. Wang et al. [18, 19] considered asymptomatic infectious individuals and contaminated environments contribute substantially to new HFMD infections. Chadsuthi and Wichapeng [20] had investigated HFMD transmission dynamics in Bangkok, Thailand, and concluded that the direct transmission from asymptomatic individuals and indirect transmission via free-living viruses are important factors to new HFMD infections.

Particularly, Zhu et al. [5] established an SEIQRS epidemic model with periodic transmission rate to investigate the spread of seasonal HFMD in Wenzhou and found that HFMD becomes an endemic disease in Wenzhou, and for controlling the spread of HFMD, it is beneficial to increase the quarantined rate or decrease the treatment cycle. In [6], the authors presented a spatial-temporal ARMA model and found that HFMD had positive spatial autocorrelation and the incidence seasonal peak was between May and July. Zhou et al. [21] investigated the epidemiological feature of HFMD in Wenzhou. Dai et al. [22] had shown that the school opening and meteorological factors were primarily responsible for annual multiple-peak pattern of the outbreaks of HFMD in Wenzhou.

But, to our knowledge, the research on the effects of the EV-A71 vaccination on the spread of HFMD in Wenzhou seems rare.

Thanks to the insightful work of Li et al. [16, 17], in this paper, we will focus on the optimal control strategies of HFMD in Wenzhou incorporating the EV-A71 vaccination. The rest of this article is organized as follows: In Section 2, we introduce the data of HFMD in Wenzhou and build a HFMD model with periodic transmission rate. In Section 3, we present the simulation results of the monitor data of Wenzhou from January 2012 to April 2018 and give the sensitivity analysis of the basic reproduction number and the optimal control strategies. In Section 4, we conclude the epidemiological significance of the results.

2. Materials and Dynamic Model

2.1. Data Source

Wenzhou is a prefecture-level city in southeastern Zhejiang province in the People’s Republic of China. It is surrounded by Yandang mountains, the East China Sea, and 436 islands, while its lowlands are almost entirely along its East China Sea coast, which is nearly 355 kilometres (221 miles) long. Most of Wenzhou’s area is mountainous as almost 76 percent of its 11,784 square kilometre (4,550 square miles) surface area is classified as mountains and hills. The 2010 China National Census registered 9.122 million residents in the Wenzhou area [21]. The prefecture-level city of Wenzhou currently administers four districts (Lucheng, Longwan, Ouhai, and Dongtou), two county-level cities (Rui’an and Yueqing), and five counties (Yongjia, Wencheng, Pingyang, Taishun, and Cangnan) (see Figure 1).

Since Wenzhou has a humid subtropical climate zone with an annual average 18.08°C (or 64.5°F), it is of particular public health significance to update molecular epidemiology of HFMD in Wenzhou [5]. HFMD data are obtained mainly from epidemiologic bulletins published by the Wenzhou Center for Disease Control and Prevention from January 2012 to April 2018 (76 months), including basic demographic characteristics of HFMD cases and daily incident cases [23]. Because Dongtou District is composed of 168 islands and the number of HFMD patients is relatively small, we only simulate the data of HFMD in Wenzhou and 10 districts’ county-level cities or counties under its jurisdiction, except Dongtou District. According to the obtained records, there are 197,821 HFMD cases.

2.2. Dynamic Model of HFMD

Suppose that the whole population is divided into six compartments at time t: susceptible with vaccination , susceptible without vaccination , exposed , clinical infectious , subclinical infectious , and recovered . Obviously, . Similar to the models in [5, 16, 17, 22], we can establish the following flowchart of compartments of HFMD (see Figure 2):

In addition, the corresponding transmission model of HFMD in Wenzhou is as follows:

All parameters in model (1) are positive, and the interpretations of them are as follows:(i)μ: natural death rate(ii)α: rate of progression to the infectious(iii): recovery rates of the infectious(iv): immune loss rates(v): recruitment rates(vi)ρ: rate of HFMD inpatients(vii)σ: the reduction in risk of infection due to EV71(viii): rates of disease transmission

It is worthy to note that, in model (1), the periodic transmission rates between , , I, and L are employed. This is indeed a well-established way of introducing seasonality into epidemic models and applied widely [2427]. The transmission rate between , , and I is taken asand the transmission rate between , , and L is taken aswhere , and are positive constants, and are the baseline contact rates, and and are the magnitudes of forcing, which can be determined by the minimum chi-square fitting [16, 28] in Section 2.3.

Since the EV-A71 vaccination was given in Wenzhou from October 2016 (the 59th month), we assume that the coefficient of reduction in exposure due to vaccination is and a piecewise function is employed to represent σ as the following:

It is easy to obtain that model (1) always has a disease-free equilibrium , where and . From model (1), we have

Therefore, the biologically feasible region for model (1) is

Easy to prove that region is positively invariant with respect to system (1).

The basic reproduction number , defined as the average number of secondary infections generated by a single infected individual introduced into a completely susceptible population [2933], is one of the important quantities in epidemiology. For model (1), following the method in [30], we can obtain as follows:where

Similar to that in [5], we can obtain the threshold dynamics of model (1) as follows.

Theorem 1. If , the disease-free equilibrium of model (1) is globally asymptotically stable, while if , model (1) has at least one positive periodic solution which is uniformly persistent and is unstable.

The proof of Theorem 1 is similar to that of Theorems 4 and 7 in [5], and we omit it here.

2.3. Parameter Estimation and Numerical Simulations

In this section, we will estimate all the parameters and initial conditions of model (1). On the basis of the biological significance of parameters, we set the upper and lower bounds of each parameter. According to the data of the sixth census, the average life expectancy of the Chinese people is 76.34 years in 2015 [34]; hence, we take the natural death rate as . According to the statistics data released by Wenzhou Municipal Bureau of Statistics [35], from 2012 to 2018, the permanent population of Wenzhou is about 9 million. Specifically, the most populous cities, Rui’an and Yueqing are about 1.4 million people every year and the cities with the smallest population (Wencheng and Taishun) are about 0.2 million people every year.

From statistics data of HFMD in Wenzhou, we can only obtain the initial value of directly (see for details Tables 1 and 2). Obviously, we cannot obtain other initial values of , , , , and . In this situation, we can take the whole population as the susceptible and assume that the upper and lower limits of the two types of the susceptible ( and ) as and for all regions and we can obtain the values of and with the help of MATLAB by using the following “multiple starting points” optimization algorithm. Also, the initial values of , , and can also be obtained in MATLAB, whose results are shown in Tables 1 and 2.

Next, we will describe the specific simulation method. Hence, we need to estimate the other 16 parameters and 5 initial values through calculating the minimum sum of chi-square [16, 28]:where shows the real value each month and shows the estimated value each month.

In order to find the optimal parameters, we choose the fminsearch function in the optimization toolbox of the software MATLAB and we try to avoid the occurrence of local optimal solutions by using the following “multiple starting points” optimization algorithm [36, 37]:(i)Step 1. Parametric hypothesis: set a reasonable range of each parameter and initial condition of model (1) to be estimated;(ii)Step 2. Evaluating simulation results: random values in the range of the first step are taken as the initial values of the first simulation. Judging whether the simulation results converge or not, the objective function value (sum of chi-square) is calculated;(iii)Step 3. Repetitive simulation: if the simulation results are not convergent, the sum of chi-square values is quite different from that of the previous simulation. Take the simulation result as a new starting value and simulate it again;(iv)Step 4. Getting a set of optimal solutions: repeat Steps 2 and 3 until the simulation results converge, and the sum of chi-square values does not change significantly;(v)Step 5. Getting multiple sets of optimal values: randomly change the starting value of the simulation and repeat Steps 2, 3, and 4. Then, we can get multiple sets of optimal parameters;(vi)Step 6. Hypothesis test: in order to screen out the results which conform to the hypothesis test, we can give the original hypothesis and alternative hypothesis, which can be seen in literatures [16, 28], and select the significant level as 0.05, . For this reason, if the chi-square value is less than this value (see Tables 1 and 2), it is considered that the original hypothesis is not rejected;(vii)Step 7. Screening out the most reasonable and optimal parameters: through the sixth step, we get several sets of optimal parameter values, further calculate the basic reproduction number, discuss the actual meaning of each parameter, and remove the unreasonable values.

Using the algorithm above, we can estimate the parameters and initial values of model (1) (see Tables 1 and 2), and the medians and arithmetic means of them are shown in Table 3 in Appendix. The range of the latent period is 1.0 d to 6.2 d, and the median value of it is about 1.7 d. The range of the course of treatment for clinical infectious is about d, and the median value is about 5.4 d. The range of the course of treatment for clinical infectious is about d, and the median value is about 2.8 d. With these parameters, we can numerically solve model (1) to fit the monthly reported HFMD cases, and the numerical results of the relations between the fitted data with the real data of Wenzhou are shown in Figure 3.

In addition, Li et al. [6] found that the distribution of HFMD in Wenzhou is heterogeneous at county level and Ouhai municipal district is the most severe region, followed by Lucheng, Longwan, and Ruian. For the sake of learning of HFMD dynamics in Wenzhou, we show the numerical results of 10 districts’ county-level cities or counties in Figure 4.

The simulation results (see Figures 3 and 4 and Tables 13) can basically describe the prevalence of HFMD in Wenzhou. Considering the vaccination again, it is true that the epidemic situation can be simulated to decline after October 2016. The pattern of HFMD dynamics is similar among the 10 districts’ county-level cities. Two peaks appear in Spring and Autumn each year in all 10 districts (see Figure 4), which is also consistent with the case of Wenzhou (see Figure 3). The estimation of is little bit larger than 1 (see Tables 1 and 2). Also, is calculated by the integral average of periodic function (see Tables 1 and 2). In Wenzhou city, . In addition, from Table 3, we can find that, except and , the medians of other parameters are similar to the arithmetic mean, which means that the incidence of HFMD in Wenzhou is basically the same as the epidemic regularity.

2.4. Sensitivity Analysis

From Theorem 1, we know that the basic reproduction number can be used to govern whether HFMD goes into extinction or not. It should be noted that there are 11 undetermined parameters in , each with different mean value and variance (see Table 3 in Appendix), which means that it is tedious and extremely complex to study the effects of these parameters on the HFMD dynamics. In order to detect the underlying factors more simply, a better idea is to do the sensitivity analysis of different categories of parameters [22].

Parametric sensitivity analysis is generally divided into two types: deterministic sensitivity analysis and uncertainty sensitivity analysis [38]. The former generally calculates the partial derivatives of the basic reproduction number with respect to each parameter, while the latter generally assumes that each parameter conforms to a certain prior distribution (e.g., uniform and standard normal distribution).

In order to identify the sensitivity of different parameters to , we used Latin hypercube sampling (LHS) and partial rank correlation coefficient (PRCC) [39] to detect the influence of each uncertain parameter on . The sample size is chosen as . All 11 parameters in the expression of are used as input values, as the output value, and the significance level is taken as 0.01. Suppose that all the input parameters are standard normal distributions, and the estimated expectations are shown in Tables 1 and 2.

The numerical results of the partial rank correlation coefficient of are shown in Table 4, and its bar chart is shown in Figure 5. The larger the PRCC in absolute value, the more important the parameters in responding to the change of new HFMD cases. Hence, we can conclude that parameters α, ρ, , , , and have positive impacts on ; conversely, , , , , and have negative impacts. Also, the most sensitive parameter for is ρ, followed by , , , , , and . These may provide potential guidance for HFMD prevention and control in strategy formulation.

3. Optimal Control Strategies

Optimal control theory deals with the problem of finding a control law for a given system such that a certain optimality criterion is achieved [16, 40, 41]. Based on the results of sensitivity analysis, we choose , , , and as controlled parameters. Considering the actual situation of HFMD epidemic in Wenzhou and the approximate medical expenditure of Wenzhou residents, we design the optimal control strategies:(i)Decreasing the transmission rates: assume that the associated force of transmission rates, and , are reduced by and , respectively. Also, and measure the precaution effort such as window ventilation, attention to disinfection, regular exercise, and HFMD prevention and control knowledge education.(ii)Increasing the recovery rates: assume that the recovery rates ( and ) are increased by and , respectively. Also, and measure the efficiency of treatment to achieve early recovery for the clinical infectious and the efficiency increasing children’s immunity (i.e., more fruits, vegetables, and regular exercise) for the subclinical infectious, respectively.

Taking into account the assumptions above, model (1) incorporating four control measures is governed by the following differential equations:with appropriate initial conditions.

For the optimal control problem of system (10), we consider the control variables relative to the state variables , , E, I, L, and R where control variables are bounded and measured withwhere T represents a given control period. Define the objective function:where represent the weight constants of the clinical infectious and subclinical infectious individuals, respectively, are weight constants for transmission rates, and and are recovery rates of infectious and , respectively.

The main focus is to minimize the number of the exposed, the clinical infectious and the subclinical infectious, and the costs required to control HFMD by using possible minimal control variables . Hence, we need to find such thatsubject to system (10).

In order to find an optimal solution, we first present the Lagrangian and Hamiltonian for the optimal control problems. Also, the Lagrangian of the problem is

We need to find the minimal value of the Lagrangian, so we define the Hamiltonian for the control problem aswhere are the adjoint functions to be determined suitably. We apply Pontryagin’s maximum principle [42]. If is an optimal solution of the optimal control problem, then there exists a nontrivial vector function satisfying the following inequalities:

Now, we apply the necessary conditions to the Hamiltonian .

Theorem 2. For problems (10)–(13) with fixed initial conditions , , , , , and and a fixed final time T, there exist adjoint functions such thatwherewith transversality conditions (boundary conditions)Furthermore, are represented by

Proof. To determine the adjoint equations and the transversality conditions, we use Hamiltonian (15). Setting , , , , , and and differentiating Hamiltonian (15) with respect to , and R, respectively, we can obtain (19). Solving the equations , on the interior of the control set and using the optimality conditions and the property of the control space U, we can derive (22).

Theorem 3. Problems (10)–(13) with given initial conditions , , , , , and are fixed final time T, which admits a unique optimal solution associated to an optimal control set on .

Proof. The existence of an optimal solution associated with an optimal control comes from the convexity of the integrand of the cost functional J with respect to the controls and the Lipschitz property of the state system with respect to state variables , and R (see [43, 44]). For small final time T, the optimal control is given by (22) that is unique by the analysis above. Because the problems (10)–(13) are autonomous, uniqueness is valid for any time T and not only for small time T.

The optimal control set predicted by Theorems 2 and 3 represents the optimal intervention strategy, given the cost constraints, and can be found by the application of celebrated Pontryagin’s maximum principle [42] and appropriate numerical methods [45].

To find out the optimal control and the state system numerically, we use the approximate algorithm for obtaining the optimal control based on the forward-backward sweep scheme which is proposed in [46]. In order to speed up the convergence, we use the convex combination as follows.

Update by entering the new adjoint variables into formula (22). They are not stored as the control variables , but as temporary vectors . The control variables are set as the convex combination of the last iteration of , namely, , and the temporary vectors . That is,where k is the current iteration and . Here, we adopt .

The control period T is 76 months, and the initial condition and parameters come from Wenzhou city in Table 2. The weight constant values in the objective functional are . The numerical settings of indicate that we hope the number of the clinical infectious will decrease more than the subclinical infectious. Also, the specific cost of each control measure is represented by the numerical values of .

We next consider three control schemes:(i)Strategy 1. Only considering to reduce contact rates and , the epidemic situation will be controlled in the 67th month (see Figure 6), which shows that hand washing and disinfection and so on can control HFMD;(ii)Strategy 2. Only considering to increase the recovery rate and , the epidemic situation will be controlled in the 22nd month (see Figure 7), which indicates that active treatment will lead to faster recovery;(iii)Strategy 3. Considering both to reduce the contact rate and and increase the recovery rate and , HFMD will be controlled in the 10th month (see Figure 8), which indicates that the optimal control strategy proposed by model (10) is effective. The time-varying optimal control parameters , , , and of Strategy 3 are shown in Figure 9.

4. Concluding Remarks

In this paper, we investigate the dynamics of a modified HFMD model (1) with the periodic transmission rate in Wenzhou, China. The value of this study lies in two aspects: mathematically, we define the reproduction number and show that it can be used to govern the stochastic dynamics of the HFMD model (1): if , the system has a unique stable disease-free equilibrium which means the extinction of HFMD; if , it has an endemic equilibrium which leads to the persistence of HFMD (see Theorem 1). In addition, we establish the optimal control strategies (see Theorem 3) of system (10) corresponding to model (1). Epidemiologically, we show that the cost of EV-A71 vaccination affects the dynamics of the model in the following aspects:(i)The EV-A71 vaccination cannot lead to the extinction of HFMD in Wenzhou: from the fitting charts 3 and 4, the use of vaccination indeed slightly decline in HFMD patients, but from the sensitivity analysis of the basic reproduction number , the cost of vaccination cannot lead to the extinction of HFMD and HFMD will periodically break out in Wenzhou. In fact, the vaccination is only useful to HFMD patients caused by EV-A71, and useless to HFMD patients caused by CX-A16 and others. In addition, according to the median of recruitment rate (see Table 4), one can obtain that , which shows that about 33.81% of newborns are vaccinated in Wenzhou. The reason of low rate of vaccination may be that EV-A71 vaccination is not included in the medical insurance plan in China. Also, the EV-A71 vaccination can reduce the incidence of HFMD slightly, so increasing the rate of EV-A71 vaccination will benefit the control of HFMD spread in Wenzhou. See Figure 10 for more details.(ii)The effect of ρ: the rate of HFMD inpatients ρ is the most sensitive parameter, which reflects that subclinical patients are the most important cause of persistent outbreaks of HFMD. This coincides with the existing literature [1517].(iii)The control strategy: we investigate the HFMD dynamics with three types of control schemes: (1) to reduce the contact rates and (see Figure 6); (2) to increase the recovery rates and (see Figure 7); (3) to reduce and and increase and (see Figure 8) and find that strategy (3) is the optimal policy. That is, in order to control the spread of HFMD, we need to pay attention to the personal hygiene to reduce the contact of the infectious persons and contaminative objects ( and ), and increasing the active treatment and immune-boosting to shorten the course of infection will help to control the epidemic of HFMD.

It should be noted that, in the numerical results in Figure 3, there are errors between the simulations and actual HFMD cases; the maximum difference between them is about 5000. In fact, model (1) is based on some assumptions and we only consider the key factors of the spread of HFMD (see Figure 2) and neglect the influences of other factors such as sudden meteorological, temperature, and humidity. Ignoring these factors can make errors of the fitting results. However, if we consider more factors in modeling, we cannot obtain any analytical results (such as Theorem 1) about model (1). On the other hand, in the present paper, we adopt two simple periodic functions (2) and (3) as the transmission rates, which are the main factor of simulation results. Also, if we choose more complex periodic functions, we will add more parameters in model (1) and we cannot obtain the main results of optimal control strategies in Theorems 2 and 3. These are the main causes to bring about the maximum difference between simulated and actual HFMD cases that is about 5000. But, from Figure 3, we can know that except several months, other months fit well and we can learn about the trend of the spread of HFMD in Wenzhou.

In future research, we will take the interference of other factors, such as randomness and intervention strategies, to obtain more practical rules of the transmission of HFMD.

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 there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (Grant nos. 61672013, 11601179, and 61772017) and Huaian Key Laboratory for Infectious Diseases Control and Prevention (HAP201704).