Optimal Control Strategies of HFMD in Wenzhou, China

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 inWenzhou, China, analytically and numerically. We define the basic reproduction number R0 and show that it can be used to determine whether HFMD becomes extinct or not. Based on the monthly reported HFMD cases inWenzhou for 76 months, we estimate the parameters in the dynamic model by using the method of minimum chisquare fitting, conduct the sensitivity analysis to investigate the influence of each uncertain parameter onR0 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 inWenzhou, we need to increase the rate of EV-A71 vaccination, decrease the contact rates, and shorten the course of disease.


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 [1][2][3]. 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].
ere 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 [9][10][11][12][13][14][15]. 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, ailand, 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.
anks 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. e 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.  [21]. e 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).

Materials and Dynamic Model
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.

Dynamic Model of HFMD.
Suppose that the whole population N(t) is divided into six compartments at time t: susceptible with vaccination S 1 (t), susceptible without vaccination S 2 (t), exposed E(t), clinical infectious I(t), subclinical infectious L(t), and recovered R(t). Obviously, N(t) � S 1 (t) + S 2 (t) + E(t) + I(t) + L(t) + R(t). 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) δ 1 and δ 2 : recovery rates of the infectious (iv) c 1 and c 2 : immune loss rates (v) Λ 1 and Λ 2 : recruitment rates (vi) ρ: rate of HFMD inpatients (vii) σ: the reduction in risk of infection due to EV71 (viii) β 1 (t) and β 2 (t): rates of disease transmission It is worthy to note that, in model (1), the periodic transmission rates between S 1 , S 2 , I, and L are employed.
is is indeed a well-established way of introducing seasonality into epidemic models and applied widely [24][25][26][27]. e transmission rate between S 1 , S 2 , and I is taken as and the transmission rate between S 1 , S 2 , and L is taken as where a 1 , b 1 , c 1 , a 2 , b 2 , and c 2 are positive constants, a 1 and a 2 are the baseline contact rates, and b 1 and b 2 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 σ 0 and a piecewise function is employed to represent σ as the following: It is easy to obtain that model (1) always has a disease- erefore, the biologically feasible region for model (1) is Easy to prove that region Ω is positively invariant with respect to system (1). e basic reproduction number R 0 , defined as the average number of secondary infections generated by a single infected individual introduced into a completely susceptible population [29][30][31][32][33], is one of the important quantities in epidemiology. For model (1), following the method in [30], we can obtain R 0 as follows: where β 1 � 1 76 Similar to that in [5], we can obtain the threshold dynamics of model (1) as follows.  Figure 1: e specific geography location of Wenzhou city, including 10 districts' county-level cities/counties [6].

Complexity
Theorem 1. If R 0 < 1, the disease-free equilibrium P 0 � (S 0 1 , S 0 2 , 0, 0, 0, 0) of model (1) is globally asymptotically stable, while if R 0 > 1, model (1) has at least one positive periodic solution which is uniformly persistent and P 0 is unstable. e proof of eorem 1 is similar to that of eorems 4 and 7 in [5], and we omit it here.

Parameter Estimation and Numerical Simulations.
In this section, we will estimate all the parameters and initial conditions of model (1) From statistics data of HFMD in Wenzhou, we can only obtain the initial value of I(0) directly (see for details Tables 1 and 2). Obviously, we cannot obtain other initial values of S 1 (0), S 2 (0), E(0), L(0), and R(0). 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 (S 1 (0) and S 2 (0)) as 10 6 and 10 4 for all regions and we can obtain the values of S 1 (0) and S 2 (0) with the help of MATLAB by using the following "multiple starting points" optimization algorithm. Also, the initial values of E(0), L(0), and R(0) 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 chisquare [16,28]: where T(t i )(i � 1, 2, . . . , 76) shows the real value each month and T(t i )(i � 1, 2, . . . , 76) 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]:  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. e range of the latent period 1/α is 1.0 d to 6.2 d, and the median value of it is about 1.7 d. e range of the course of treatment for clinical infectious 1/δ 1 is about 2.9 − 10.0 d, and the median value is about 5.4 d. e range of the course of treatment for clinical infectious 1/δ 2 is about 2.0 − 9.5 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   Complexity 5 Figure 4), which is also consistent with the case of Wenzhou (see Figure 3). e estimation of R 0 is little bit larger than 1 (see Tables 1 and 2). Also, R 0 is calculated by the integral average of periodic function (see Tables 1 and 2). In Wenzhou city, R 0 � 1.0064 > 1. In addition, from Table 3, we can find that, except Λ 1 and c 1 , 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.

Sensitivity Analysis.
From eorem 1, we know that the basic reproduction number R 0 can be used to govern whether HFMD goes into extinction or not. It should be noted that there are 11 undetermined parameters in R 0 , 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]. e 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 R 0 , we used Latin hypercube sampling (LHS) and partial rank correlation coefficient (PRCC) [39] to detect the influence of each uncertain parameter on R 0 . e sample size is chosen as n � 2000. All 11 parameters in the expression of R 0 are used as input values, R 0 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. e numerical results of the partial rank correlation coefficient of R 0 are shown in Table 4, and its bar chart is shown in Figure 5. e 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 α, ρ, β 1 , β 2 , β 1 , and β 2 have positive impacts on R 0 ; conversely, δ 1 , δ 2 , c 1 , c 2 , and σ 0 have negative impacts. Also, the most sensitive parameter for R 0 is ρ, followed by δ 1 , β 2 , β 1 , β 1 , δ 2 , and β 2 . ese may provide potential guidance for HFMD prevention and control in strategy formulation.

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 β 1 (t), β 2 (t), δ 1 , and δ 2 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, β 1 (t) and β 2 (t), are reduced by u 1 (t) and u 2 (t), respectively. Also, 1 − u 1 (t) and 1 − u 2 (t) 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 (δ 1 and δ 2 ) are increased by u 3 (t) and u 4 (t), respectively. Also, 1 + u 3 (t) and 1 + u 4 (t) measure the efficiency of treatment to achieve early recovery for the clinical infectious I(t) 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:    Complexity 8 Complexity with appropriate initial conditions. For the optimal control problem of system (10), we consider the control variables u(t) � (u 1 , u 2 , u 3 , u 4 ) ∈ U relative to the state variables S 1 , S 2 , E, I, L, and R where control variables are bounded and measured with where T represents a given control period. Define the objective function: where A i (i � 1, 2) represent the weight constants of the clinical infectious and subclinical infectious individuals, respectively, B i (i � 1, 2, 3, 4) are weight constants for transmission rates, and β 1 (t) and β 2 (t) are recovery rates of infectious δ 1 and δ 2 , respectively. e 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 u i (i � 1, 2, 3, 4). Hence, we need to find (u * 1 , u * 2 , u * 3 , u * 4 ) such that subject 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 H for the control problem as where λ i (i � 1, 2, . . . , 6) are the adjoint functions to be determined suitably. We apply Pontryagin's maximum principle [42]. If (x, u) is an optimal solution of the optimal control problem, then there exists a nontrivial vector function λ � (λ 1 , λ 2 , λ 3 , λ 4 , λ 5 , λ 6 ) satisfying the following inequalities: Now, we apply the necessary conditions to the Hamiltonian H.

Theorem 2.
For problems (10)- (13) with fixed initial conditions S 1 (0), S 2 (0), E(0), I(0), L(0), and R(0) and a fixed final time T, there exist adjoint functions λ i (i � 1, 2, . . . , 6) such that the Lipschitz property of the state system with respect to state variables S 1 , S 2 , E, I, L, 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.
□ e optimal control set predicted by eorems 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 forwardbackward sweep scheme which is proposed in [46]. In order to speed up the convergence, we use the convex combination as follows.
Update u i by entering the new adjoint variables λ i , i � 1, 2, . . . , 6, into formula (22). ey are not stored as the control variables u i , but as temporary vectors u i . e control variables u i are set as the convex combination of the last iteration of u i , namely, oldu i , and the temporary vectors u i . at is, where k is the current iteration and c ∈ (0, 1). Here, we adopt c � 0.2.
e control period T is 76 months, and the initial condition and parameters come from Wenzhou city in Table 2 We next consider three control schemes: (i) Strategy 1. Only considering to reduce contact rates β 1 (t) and β 2 (t), 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 δ 1 and δ 2 , 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 β 1 (t) and β 2 (t) and increase the recovery rate δ 1 and δ 2 , HFMD will be controlled in the 10th month (see Figure 8), which indicates that the optimal control strategy proposed by model (10) is effective. e timevarying optimal control parameters u 1 , u 2 , u 3 , and u 4 of Strategy 3 are shown in Figure 9.

Concluding Remarks
In this paper, we investigate the dynamics of a modified HFMD model (1) with the periodic transmission rate in Wenzhou, China. e value of this study lies in two aspects:   mathematically, we define the reproduction number R 0 and show that it can be used to govern the stochastic dynamics of the HFMD model (1): if R 0 < 1, the system has a unique stable disease-free equilibrium which means the extinction of HFMD; if R 0 > 1, it has an endemic equilibrium which leads to the persistence of HFMD (see eorem 1). In addition, we establish the optimal control strategies (see eorem 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) e 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 R 0 , 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 Λ 1 /Λ 1 + Λ 2 � 33.81%, which shows that about 33.81% of newborns are vaccinated in Wenzhou. e 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) e 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. is coincides with the existing literature [15][16][17]. (iii) e control strategy: we investigate the HFMD dynamics with three types of control schemes: (1) to reduce the contact rates β 1 (t) and β 2 (t) (see Figure 6); (2) to increase the recovery rates δ 1 and δ 2 (see Figure 7); (3) to reduce β 1 (t) and β 2 (t) and increase δ 1 and δ 2 (see Figure 8) and find that strategy (3) is the optimal policy. at 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 (β 1 and β 2 ), 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 eorem 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 eorems 2 and 3. ese 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
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.