# Optimal Control Strategies of HFMD in Wenzhou, China

Authors: Zuqin Ding; Yong Li; Yongli Cai; Yueping Dong; Weiming Wang
Journal: Complexity
 (2020)

---

## 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 numberR0 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 R0 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.


---

## Body

## 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 [1–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].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 [9–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, 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).Figure 1
The specific geography location of Wenzhou city, including 10 districts’ county-level cities/counties [6].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 populationNt is divided into six compartments at time t: susceptible with vaccination S1t, susceptible without vaccination S2t, exposed Et, clinical infectious It, subclinical infectious Lt, and recovered Rt. Obviously, Nt=S1t+S2t+Et+It+Lt+Rt. Similar to the models in [5, 16, 17, 22], we can establish the following flowchart of compartments of HFMD (see Figure 2):Figure 2
Flowchart of compartments of the HFMD model.In addition, the corresponding transmission model of HFMD in Wenzhou is as follows:(1)dS1dt=Λ1−σβ1tS1IN−σβ2tS1LN+γ1R−μS1,dS2dt=Λ2−β1tS2IN−β2tS2LN+γ2R−μS2,dEdt=σβ1tS1IN+σβ2tS1LN+β1tS2IN+β2tS2LN−α+μE,dIdt=ραE−δ1+μI,dLdt=1−ραE−δ2+μL,dRdt=δ1I+δ2L−γ1+γ2+μR.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)
δ1andδ2: recovery rates of the infectious(iv)
γ1andγ2: immune loss rates(v)
Λ1andΛ2: recruitment rates(vi)
ρ: rate of HFMD inpatients(vii)
σ: the reduction in risk of infection due to EV71(viii)
β1tandβ2t: rates of disease transmissionIt is worthy to note that, in model (1), the periodic transmission rates between S1, S2, I, and L are employed. This is indeed a well-established way of introducing seasonality into epidemic models and applied widely [24–27]. The transmission rate between S1, S2, and I is taken as(2)β1t≔a1+1+b1sinπtT+c1,and the transmission rate between S1, S2, and L is taken as(3)β2t≔a2+1+b2sinπtT+c2,where a1,b1,c1,a2,b2, and c2 are positive constants, a1 and a2 are the baseline contact rates, and b1 and b2 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:(4)σ=1,1≤t≤58,σ0,59≤t≤76.It is easy to obtain that model (1) always has a disease-free equilibrium P0=S10,S20,0,0,0,0, where S10=Λ1/μ and S20=Λ2/μ. From model (1), we have(5)dNdt=Λ1+Λ2−μN.Therefore, the biologically feasible region for model (1) is(6)Ω=S1,S2,E,I,L,R∈R+60<S1+S2+E+I+L+R≤Λ1+Λ2μ.Easy to prove that regionΩ is positively invariant with respect to system (1).The basic reproduction numberR0, defined as the average number of secondary infections generated by a single infected individual introduced into a completely susceptible population [29–33], is one of the important quantities in epidemiology. For model (1), following the method in [30], we can obtain R0 as follows:(7)R0≔β˜1S10S10+S20+β^1S20S10+S20ραα+μδ1+μ+β˜2S10S10+S20+β^2S20S10+S201−ραα+μδ2+μ,where(8)β˜1=176∫158β1tdt+σ0∫5976β1tdt,β˜2=176∫158β2tdt+σ0∫5976β2tdt,β^1=176∫176β1tdt,β^2=176∫176β2tdt.Similar to that in [5], we can obtain the threshold dynamics of model (1) as follows.Theorem 1.
IfR0<1, the disease-free equilibrium P0=S10,S20,0,0,0,0 of model (1) is globally asymptotically stable, while if R0>1, model (1) has at least one positive periodic solution which is uniformly persistent and P0 is unstable.The proof of Theorem1 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 μ=1/76.34×12=1.0916×10−2. 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 ofI0 directly (see for details Tables 1 and 2). Obviously, we cannot obtain other initial values of S10, S20, E0, L0, and R0. 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 (S10 and S20) as 106 and 104 for all regions and we can obtain the values of S10 and S20 with the help of MATLAB by using the following “multiple starting points” optimization algorithm. Also, the initial values of E0, L0, and R0 can also be obtained in MATLAB, whose results are shown in Tables 1 and 2.Table 1
Parameters and their values.
Range of parameter
Wenzhou
Lucheng
Longwan
Ouhai
Rui’an
α
∈
0,30
30.0000
14.7784
20.6481
6.1205
6.2197
δ
1
∈
3,15
4.4101
5.4885
3.8918
5.7724
10.3272
δ
2
∈
3,15
7.0402
12.1971
11.2610
11.0008
7.0558
γ
1
∈
0,10
0.0112
5.2515
0.2798
3.3819
0.2913
γ
2
∈
0,10
4.2300
9.8776
9.9723
9.5839
9.9936
Λ
1
∈
0,100
0.0001
29.5983
8.01
×
10
−
5
16.7320
6.31
×
10
−
6
Λ
2
∈
0,100
0.0002
29.3058
0.0106
7.6879
32.9412
ρ
∈
0,0.5
0.1610
0.1670
0.1346
0.1747
0.3090
σ
0
∈
0,0.5
0.7473
0.9262
0.1146
0.9438
0.8278
T
∈
6,18
6.1037
6.1231
6.0470
6.1428
6.0046
a
1
∈
0,30
23.897
29.9995
25.6717
30.0000
29.5449
a
2
∈
0,30
6.17
×
10
−
8
0.1291
0.0753
0.1256
0.0002
b
1
∈
−
1,1
0.9604
0.7050
0.9110
0.9135
0.5348
b
2
∈
−
1,1
0.7911
1.0000
0.9997
1.0000
0.7523
c
1
∈
0,1
0.0532
0.1160
0.0012
0.1665
0.0670
c
2
∈
0,1
0.2992
0.1284
0.1507
0.0464
1.09
×
10
−
6
β
˜
1
23.1928
30.1031
21.0575
30.2313
28.9888
β
˜
2
0.9568
1.1290
0.8659
1.1262
0.9792
β
^
1
24.9188
31.0146
26.6997
31.0177
30.5655
β
^
2
1.0223
1.1508
1.1093
1.1423
1.0281
S
1
0
∈
10
4
,
10
6
863190
130216
83540
109356
335684
S
2
0
∈
10
4
,
10
6
995989
107109
999996
119795
10
6
E
0
∈
0
,
10
5
13
564
104
1378
0
I
0
1091
89
38
39
215
L
0
∈
0
,
10
5
1
0
0
3319
12051
R
0
∈
0
,
10
6
458
6818
10241
4259
0
χ
2
—
35.3842
0.8659
25.2958
17.7474
R
0
1.0015
1.0042
1.0034
1.0030
1.0116Table 2
Parameters and their values.
Range of parameters
Yongjia
Yueqing
Wencheng
Pingyang
Taishun
Cangnan
α
∈
0,30
4.9592
29.9046
27.7359
27.8485
29.0212
15.0574
δ
1
∈
3,15
8.1110
4.3826
4.2318
5.3634
3.0259
5.9019
δ
2
∈
3,15
3.2000
14.9895
13.2424
12.2123
14.4614
7.9358
γ
1
∈
0,10
3.85
×
10
−
5
0.2154
2.83
×
10
−
6
2.0229
0.0225
5.34
×
10
−
6
γ
2
∈
0,10
7.4057
9.9999
8.3133
8.6110
8.3528
9.9993
Λ
1
∈
0,100
0.0006
6.82
×
10
−
5
2.35
×
10
−
5
9.6393
3.40
×
10
−
6
69.2267
Λ
2
∈
0,100
0.0217
0.0524
41.3280
99.9888
0.0513
33.7068
ρ
∈
0,0.5
0.2402
0.1970
0.1599
0.1665
0.1222
0.1718
σ
0
∈
0,0.5
0.5428
0.4965
0.9104
0.9950
0.8878
0.9862
T
∈
6,18
6.0000
6.0755
6.0239
6.0679
6.2586
6.0357
a
1
∈
0,30
25.2706
20.1875
23.6496
29.3384
22.3592
29.9524
a
2
∈
0,30
0.0004
0.1394
0.1961
0.1180
0.0064
0.1250
b
1
∈
−
1,1
0.6885
0.9267
0.4246
0.9999
0.5025
0.9516
b
2
∈
−
1,1
0.4533
0.9998
0.6541
0.5145
0.8819
0.8654
c
1
∈
0,1
0.1425
0.0667
0.1098
0.0940
0.1448
6.86
×
10
−
7
c
2
∈
0,1
0.2421
0.1338
0.1972
0.1625
0.0776
0.1284
β
˜
1
23.2590
18.5501
23.8494
29.9458
22.4753
30.4902
β
˜
2
0.8983
1.0254
1.1848
1.1235
0.9842
1.1453
β
^
1
26.2984
21.2133
24.6651
30.3681
23.3613
30.9833
β
^
2
1.0194
1.1687
1.2211
1.1341
1.0095
1.1556
S
1
0
∈
10
4
,
10
6
25732
205161
144028
78907
269882
99081
S
2
0
∈
10
4
,
10
6
971958
223628
183048
136866
10
6
108336
E
0
∈
0
,
10
5
28
148
117
925
89
831
I
0
93
228
19
204
18
132
L
0
∈
0
,
10
5
3450
0
1
1
0
860
R
0
∈
0
,
10
6
1875
7037
5352
26
7589
3708
χ
2
4.8267
49.6591
1.1284
31.9636
0.0831
22.7046
R
0
1.0124
1.0128
1.0059
1.0163
1.0004
1.0089Next, 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]:(9)Jθ=∑i=176Tti−T^ti2T^ti,where Ttii=1,2,…,76 shows the real value each month and T^tii=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]:(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, χ0.05254=72.153. 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 1/α 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 1/δ1 is about 2.9−10.0 d, and the median value is about 5.4 d. The 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.Table 3
The medians and arithmetic means of the parameters.
Parameter
Median
Arithmetic mean
Source
α
∈
0,30
18.2294
25.3815
Estimate
δ
1
∈
3,15
5.6497
5.4260
Estimate
δ
2
∈
3,15
10.7556
12.2047
Estimate
γ
1
∈
0,10
1.1465
0.2856
Estimate
γ
2
∈
0,10
9.2109
9.7308
Estimate
Λ
1
∈
0,100
12.5196
0.0003
Estimate
Λ
2
∈
0,100
24.5073
18.4969
Estimate
ρ
∈
0,0.5
0.1843
0.1733
Estimate
σ
0
∈
0,0.5
0.7631
0.8991
Estimate
T
∈
6,18
6.0779
6.0575
Estimate
β
˜
1
25.8951
26.4191
Calculate
β
˜
2
1.0462
1.0745
Calculate
β
^
1
27.6188
28.5339
Calculate
β
^
2
1.1139
1.1382
Calculate
R
0
1.0079
1.0103
CalculateFigure 3
The reported data vs. the numerical solution of the HFMD model (1) in Wenzhou.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.Figure 4
The reported data vs. the numerical solution of the HFMD model (1) in (a) Lucheng, (b) Ouhai, (c) Longwan, (d) Rui’an, (e) Cangnan, (f) Pingyang, (g) Yueqing, (h) Yongjia, (i) Taishun, and (j) Wencheng.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)The simulation results (see Figures3 and 4 and Tables 1–3) 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 R0 is little bit larger than 1 (see Tables 1 and 2). Also, R0 is calculated by the integral average of periodic function (see Tables 1 and 2). In Wenzhou city, R0˜=1.0064>1. In addition, from Table 3, we can find that, except Λ1 and γ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.


### 2.4. Sensitivity Analysis

From Theorem1, we know that the basic reproduction number R0 can be used to govern whether HFMD goes into extinction or not. It should be noted that there are 11 undetermined parameters in R0, 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 toR0, we used Latin hypercube sampling (LHS) and partial rank correlation coefficient (PRCC) [39] to detect the influence of each uncertain parameter on R0. The sample size is chosen as n=2000. All 11 parameters in the expression of R0 are used as input values, R0 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 ofR0 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 α, ρ, β˜1, β˜2, β^1, and β^2 have positive impacts on R0; conversely, δ1, δ2, γ1, γ2, and σ0 have negative impacts. Also, the most sensitive parameter for R0 is ρ, followed by δ1, β^2, β^1, β˜1, δ2, and β˜2. These may provide potential guidance for HFMD prevention and control in strategy formulation.Table 4
PRCC values forR0.
Input parameter
PRCC
P value
α
0.039809769392456
0.075818905660133
δ
1
−
0.845417442920200
0
δ
2
−
0.233470793181710
0
γ
1
−
0.013915647248970
0.534988619354773
γ
2
−
0.013936543241846
0.534375545613689
ρ
0.992295606945784
0
σ
0
−
0.036960654931174
0.099286931216591
β
˜
1
0.468625053788954
0
β
˜
2
0.169141046318725
0.000000000000031
β
^
1
0.657754383903920
0
β
^
2
0.677554835647717
0Figure 5
The values of PRCC on the outcome ofR0. All parameter values are derived from Wenzhou city in Table 2.


## 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).Figure 1
The specific geography location of Wenzhou city, including 10 districts’ county-level cities/counties [6].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 populationNt is divided into six compartments at time t: susceptible with vaccination S1t, susceptible without vaccination S2t, exposed Et, clinical infectious It, subclinical infectious Lt, and recovered Rt. Obviously, Nt=S1t+S2t+Et+It+Lt+Rt. Similar to the models in [5, 16, 17, 22], we can establish the following flowchart of compartments of HFMD (see Figure 2):Figure 2
Flowchart of compartments of the HFMD model.In addition, the corresponding transmission model of HFMD in Wenzhou is as follows:(1)dS1dt=Λ1−σβ1tS1IN−σβ2tS1LN+γ1R−μS1,dS2dt=Λ2−β1tS2IN−β2tS2LN+γ2R−μS2,dEdt=σβ1tS1IN+σβ2tS1LN+β1tS2IN+β2tS2LN−α+μE,dIdt=ραE−δ1+μI,dLdt=1−ραE−δ2+μL,dRdt=δ1I+δ2L−γ1+γ2+μR.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)
δ1andδ2: recovery rates of the infectious(iv)
γ1andγ2: immune loss rates(v)
Λ1andΛ2: recruitment rates(vi)
ρ: rate of HFMD inpatients(vii)
σ: the reduction in risk of infection due to EV71(viii)
β1tandβ2t: rates of disease transmissionIt is worthy to note that, in model (1), the periodic transmission rates between S1, S2, I, and L are employed. This is indeed a well-established way of introducing seasonality into epidemic models and applied widely [24–27]. The transmission rate between S1, S2, and I is taken as(2)β1t≔a1+1+b1sinπtT+c1,and the transmission rate between S1, S2, and L is taken as(3)β2t≔a2+1+b2sinπtT+c2,where a1,b1,c1,a2,b2, and c2 are positive constants, a1 and a2 are the baseline contact rates, and b1 and b2 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:(4)σ=1,1≤t≤58,σ0,59≤t≤76.It is easy to obtain that model (1) always has a disease-free equilibrium P0=S10,S20,0,0,0,0, where S10=Λ1/μ and S20=Λ2/μ. From model (1), we have(5)dNdt=Λ1+Λ2−μN.Therefore, the biologically feasible region for model (1) is(6)Ω=S1,S2,E,I,L,R∈R+60<S1+S2+E+I+L+R≤Λ1+Λ2μ.Easy to prove that regionΩ is positively invariant with respect to system (1).The basic reproduction numberR0, defined as the average number of secondary infections generated by a single infected individual introduced into a completely susceptible population [29–33], is one of the important quantities in epidemiology. For model (1), following the method in [30], we can obtain R0 as follows:(7)R0≔β˜1S10S10+S20+β^1S20S10+S20ραα+μδ1+μ+β˜2S10S10+S20+β^2S20S10+S201−ραα+μδ2+μ,where(8)β˜1=176∫158β1tdt+σ0∫5976β1tdt,β˜2=176∫158β2tdt+σ0∫5976β2tdt,β^1=176∫176β1tdt,β^2=176∫176β2tdt.Similar to that in [5], we can obtain the threshold dynamics of model (1) as follows.Theorem 1.
IfR0<1, the disease-free equilibrium P0=S10,S20,0,0,0,0 of model (1) is globally asymptotically stable, while if R0>1, model (1) has at least one positive periodic solution which is uniformly persistent and P0 is unstable.The proof of Theorem1 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 μ=1/76.34×12=1.0916×10−2. 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 ofI0 directly (see for details Tables 1 and 2). Obviously, we cannot obtain other initial values of S10, S20, E0, L0, and R0. 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 (S10 and S20) as 106 and 104 for all regions and we can obtain the values of S10 and S20 with the help of MATLAB by using the following “multiple starting points” optimization algorithm. Also, the initial values of E0, L0, and R0 can also be obtained in MATLAB, whose results are shown in Tables 1 and 2.Table 1
Parameters and their values.
Range of parameter
Wenzhou
Lucheng
Longwan
Ouhai
Rui’an
α
∈
0,30
30.0000
14.7784
20.6481
6.1205
6.2197
δ
1
∈
3,15
4.4101
5.4885
3.8918
5.7724
10.3272
δ
2
∈
3,15
7.0402
12.1971
11.2610
11.0008
7.0558
γ
1
∈
0,10
0.0112
5.2515
0.2798
3.3819
0.2913
γ
2
∈
0,10
4.2300
9.8776
9.9723
9.5839
9.9936
Λ
1
∈
0,100
0.0001
29.5983
8.01
×
10
−
5
16.7320
6.31
×
10
−
6
Λ
2
∈
0,100
0.0002
29.3058
0.0106
7.6879
32.9412
ρ
∈
0,0.5
0.1610
0.1670
0.1346
0.1747
0.3090
σ
0
∈
0,0.5
0.7473
0.9262
0.1146
0.9438
0.8278
T
∈
6,18
6.1037
6.1231
6.0470
6.1428
6.0046
a
1
∈
0,30
23.897
29.9995
25.6717
30.0000
29.5449
a
2
∈
0,30
6.17
×
10
−
8
0.1291
0.0753
0.1256
0.0002
b
1
∈
−
1,1
0.9604
0.7050
0.9110
0.9135
0.5348
b
2
∈
−
1,1
0.7911
1.0000
0.9997
1.0000
0.7523
c
1
∈
0,1
0.0532
0.1160
0.0012
0.1665
0.0670
c
2
∈
0,1
0.2992
0.1284
0.1507
0.0464
1.09
×
10
−
6
β
˜
1
23.1928
30.1031
21.0575
30.2313
28.9888
β
˜
2
0.9568
1.1290
0.8659
1.1262
0.9792
β
^
1
24.9188
31.0146
26.6997
31.0177
30.5655
β
^
2
1.0223
1.1508
1.1093
1.1423
1.0281
S
1
0
∈
10
4
,
10
6
863190
130216
83540
109356
335684
S
2
0
∈
10
4
,
10
6
995989
107109
999996
119795
10
6
E
0
∈
0
,
10
5
13
564
104
1378
0
I
0
1091
89
38
39
215
L
0
∈
0
,
10
5
1
0
0
3319
12051
R
0
∈
0
,
10
6
458
6818
10241
4259
0
χ
2
—
35.3842
0.8659
25.2958
17.7474
R
0
1.0015
1.0042
1.0034
1.0030
1.0116Table 2
Parameters and their values.
Range of parameters
Yongjia
Yueqing
Wencheng
Pingyang
Taishun
Cangnan
α
∈
0,30
4.9592
29.9046
27.7359
27.8485
29.0212
15.0574
δ
1
∈
3,15
8.1110
4.3826
4.2318
5.3634
3.0259
5.9019
δ
2
∈
3,15
3.2000
14.9895
13.2424
12.2123
14.4614
7.9358
γ
1
∈
0,10
3.85
×
10
−
5
0.2154
2.83
×
10
−
6
2.0229
0.0225
5.34
×
10
−
6
γ
2
∈
0,10
7.4057
9.9999
8.3133
8.6110
8.3528
9.9993
Λ
1
∈
0,100
0.0006
6.82
×
10
−
5
2.35
×
10
−
5
9.6393
3.40
×
10
−
6
69.2267
Λ
2
∈
0,100
0.0217
0.0524
41.3280
99.9888
0.0513
33.7068
ρ
∈
0,0.5
0.2402
0.1970
0.1599
0.1665
0.1222
0.1718
σ
0
∈
0,0.5
0.5428
0.4965
0.9104
0.9950
0.8878
0.9862
T
∈
6,18
6.0000
6.0755
6.0239
6.0679
6.2586
6.0357
a
1
∈
0,30
25.2706
20.1875
23.6496
29.3384
22.3592
29.9524
a
2
∈
0,30
0.0004
0.1394
0.1961
0.1180
0.0064
0.1250
b
1
∈
−
1,1
0.6885
0.9267
0.4246
0.9999
0.5025
0.9516
b
2
∈
−
1,1
0.4533
0.9998
0.6541
0.5145
0.8819
0.8654
c
1
∈
0,1
0.1425
0.0667
0.1098
0.0940
0.1448
6.86
×
10
−
7
c
2
∈
0,1
0.2421
0.1338
0.1972
0.1625
0.0776
0.1284
β
˜
1
23.2590
18.5501
23.8494
29.9458
22.4753
30.4902
β
˜
2
0.8983
1.0254
1.1848
1.1235
0.9842
1.1453
β
^
1
26.2984
21.2133
24.6651
30.3681
23.3613
30.9833
β
^
2
1.0194
1.1687
1.2211
1.1341
1.0095
1.1556
S
1
0
∈
10
4
,
10
6
25732
205161
144028
78907
269882
99081
S
2
0
∈
10
4
,
10
6
971958
223628
183048
136866
10
6
108336
E
0
∈
0
,
10
5
28
148
117
925
89
831
I
0
93
228
19
204
18
132
L
0
∈
0
,
10
5
3450
0
1
1
0
860
R
0
∈
0
,
10
6
1875
7037
5352
26
7589
3708
χ
2
4.8267
49.6591
1.1284
31.9636
0.0831
22.7046
R
0
1.0124
1.0128
1.0059
1.0163
1.0004
1.0089Next, 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]:(9)Jθ=∑i=176Tti−T^ti2T^ti,where Ttii=1,2,…,76 shows the real value each month and T^tii=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]:(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, χ0.05254=72.153. 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 1/α 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 1/δ1 is about 2.9−10.0 d, and the median value is about 5.4 d. The 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.Table 3
The medians and arithmetic means of the parameters.
Parameter
Median
Arithmetic mean
Source
α
∈
0,30
18.2294
25.3815
Estimate
δ
1
∈
3,15
5.6497
5.4260
Estimate
δ
2
∈
3,15
10.7556
12.2047
Estimate
γ
1
∈
0,10
1.1465
0.2856
Estimate
γ
2
∈
0,10
9.2109
9.7308
Estimate
Λ
1
∈
0,100
12.5196
0.0003
Estimate
Λ
2
∈
0,100
24.5073
18.4969
Estimate
ρ
∈
0,0.5
0.1843
0.1733
Estimate
σ
0
∈
0,0.5
0.7631
0.8991
Estimate
T
∈
6,18
6.0779
6.0575
Estimate
β
˜
1
25.8951
26.4191
Calculate
β
˜
2
1.0462
1.0745
Calculate
β
^
1
27.6188
28.5339
Calculate
β
^
2
1.1139
1.1382
Calculate
R
0
1.0079
1.0103
CalculateFigure 3
The reported data vs. the numerical solution of the HFMD model (1) in Wenzhou.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.Figure 4
The reported data vs. the numerical solution of the HFMD model (1) in (a) Lucheng, (b) Ouhai, (c) Longwan, (d) Rui’an, (e) Cangnan, (f) Pingyang, (g) Yueqing, (h) Yongjia, (i) Taishun, and (j) Wencheng.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)The simulation results (see Figures3 and 4 and Tables 1–3) 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 R0 is little bit larger than 1 (see Tables 1 and 2). Also, R0 is calculated by the integral average of periodic function (see Tables 1 and 2). In Wenzhou city, R0˜=1.0064>1. In addition, from Table 3, we can find that, except Λ1 and γ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.


## 2.4. Sensitivity Analysis

From Theorem1, we know that the basic reproduction number R0 can be used to govern whether HFMD goes into extinction or not. It should be noted that there are 11 undetermined parameters in R0, 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 toR0, we used Latin hypercube sampling (LHS) and partial rank correlation coefficient (PRCC) [39] to detect the influence of each uncertain parameter on R0. The sample size is chosen as n=2000. All 11 parameters in the expression of R0 are used as input values, R0 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 ofR0 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 α, ρ, β˜1, β˜2, β^1, and β^2 have positive impacts on R0; conversely, δ1, δ2, γ1, γ2, and σ0 have negative impacts. Also, the most sensitive parameter for R0 is ρ, followed by δ1, β^2, β^1, β˜1, δ2, and β˜2. These may provide potential guidance for HFMD prevention and control in strategy formulation.Table 4
PRCC values forR0.
Input parameter
PRCC
P value
α
0.039809769392456
0.075818905660133
δ
1
−
0.845417442920200
0
δ
2
−
0.233470793181710
0
γ
1
−
0.013915647248970
0.534988619354773
γ
2
−
0.013936543241846
0.534375545613689
ρ
0.992295606945784
0
σ
0
−
0.036960654931174
0.099286931216591
β
˜
1
0.468625053788954
0
β
˜
2
0.169141046318725
0.000000000000031
β
^
1
0.657754383903920
0
β
^
2
0.677554835647717
0Figure 5
The values of PRCC on the outcome ofR0. All parameter values are derived from Wenzhou city in Table 2.


## 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 β1t, β2t, δ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,β1t and β2t, are reduced by u1t and u2t, respectively. Also, 1−u1t and 1−u2t 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 u3t and u4t, respectively. Also, 1+u3t and 1+u4t measure the efficiency of treatment to achieve early recovery for the clinical infectious It 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:(10)dS1dt=Λ1−σβ1t1−u1S1IN−σβ2t1−u2S1LN+γ1R−μS1,dS2dt=Λ2−β1t1−u1S2IN−σβ2t1−u2S2LN+γ1R−μS2,dEdt=σβ1t1−u1S1IN+σβ2t1−u2S1LN+β1t1−u1S2IN+σβ2t1−u2S2LN−α+μE,dIdt=ραE−δ11+u3+μI,dLdt=1−ραE−δ21+u4+μL,dRdt=δ11+u3+μI+δ21+u4+μL−γ1+γ2+μR,with appropriate initial conditions.For the optimal control problem of system (10), we consider the control variables ut=u1,u2,u3,u4∈U relative to the state variables S1, S2, E, I, L, and R where control variables are bounded and measured with(11)U=u1,u2,u3,u4ui∈L0,1,0≤ui≤1,t∈0,T,i=1,2,3,4,where T represents a given control period. Define the objective function:(12)Ju1,u2,u3,u4=∫0TA1I+A2L+∑i=1412Biui2dt,where Aii=1,2 represent the weight constants of the clinical infectious and subclinical infectious individuals, respectively, Bii=1,2,3,4 are weight constants for transmission rates, and β1t and β2t are recovery rates of infectious δ1 and δ2, 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 variablesuii=1,2,3,4. Hence, we need to find u1∗,u2∗,u3∗,u4∗ such that(13)Ju1∗,u2∗,u3∗,u4∗=minUJu1,u2,u3,u4,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(14)LI,L,u1,u2,u3,u4=A1I+A2L+∑i=1412Biui2.We need to find the minimal value of the Lagrangian, so we define the HamiltonianH for the control problem as(15)H=LI,L,u1,u2,u3,u4+λ1dS1dt+λ2dS2dt+λ3dEdt+λ4dIdt+λ5dLdt+λ6dRdt,where λii=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:(16)dxdt=∂Ht,x,u,λ∂λ,(17)0=∂Ht,x,u,λ∂u,(18)dλdt=−∂Ht,x,u,λ∂x.Now, we apply the necessary conditions to the HamiltonianH.Theorem 2.
For problems (10)–(13) with fixed initial conditions S10, S20, E0, I0, L0, and R0 and a fixed final time T, there exist adjoint functions λii=1,2,…,6 such that(19)dλ1dt=λ1−λ3f1+f2+λ3−λ2f3+f4+λ1μ,dλ2dt=λ3−λ1f5+f6+λ2−λ3f7+f8+λ2μ,dλ3dt=λ3−λ1f5+f6+λ3−λ2f3+f4+λ3α+μ−λ4ρα−λ51−ρα,dλ4dt=−A1+λ1−λ3f9−f6+λ2−λ3f10−f4+λ4δ11+u3+μ−λ6δ11+u3,dλ5dt=−A2+λ3−λ1f5−f11+λ3−λ2f3−f12+λ5δ21+u4+μ−λ6δ21+u4,dλ6dt=λ3−λ1f5+f6+λ3−λ2f3+f4+λ6γ1+γ2+μ,where(20)f1=σβ1t1−u1IN−S1N2,f2=σβ2t1−u2LN−S1N2,f3=β1t1−u1S2IN2,f4=β2t1−u2S2LN2,f5=σβ1t1−u1S1IN2,f6=σβ2t1−u2S1LN2,f7=β1t1−u1IN−S2N2,f8=β2t1−u2LN−S2N2,f9=σβ1t1−u1S1N−IN2,f10=β1t1−u1S2N−IN2,f11=σβ2t1−u2S1N−LN2,f12=β2t1−u2S2N−LN2,with transversality conditions (boundary conditions)(21)λiT=0,i=1,2,…,6.
Furthermore,ui∗i=1,2,3,4 are represented by(22)u1∗=maxmin1,λ3−λ1σβ1tS1∗I∗+λ3−λ2β1tS2∗I∗B1N∗,0,u2∗=maxmin1,λ3−λ1σβ2tS1∗L∗+λ3−λ2β2tS2∗L∗B2N∗,0,u3∗=maxmin1,λ4−λ6δ1I∗B3,0,u4∗=maxmin1,λ5−λ6δ2L∗B4,0.Proof.
To determine the adjoint equations and the transversality conditions, we use Hamiltonian (15). Setting S1t=S1∗t, S2t=S2∗t, Et=E∗t, It=I∗t, Lt=L∗t, and Rt=R∗t and differentiating Hamiltonian (15) with respect to S1,S2,E,I,L, and R, respectively, we can obtain (19). Solving the equations ∂H/∂ui=0,i=1,2,3, 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 S10, S20, E0, I0, L0, and R0 are fixed final time T, which admits a unique optimal solution S1∗,S2∗,E∗,I∗,L∗,R∗ associated to an optimal control set u1∗,u2∗,u3∗,u4∗ on 0,T.Proof.
The existence of an optimal solutionS1∗,S2∗,E∗,I∗,L∗,R∗ associated with an optimal control u1∗,u2∗,u3∗,u4∗ comes from the convexity of the integrand of the cost functional J with respect to the controls u1,u2,u3,u4 and the Lipschitz property of the state system with respect to state variables S1,S2,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.The optimal control set predicted by Theorems2 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.Updateui by entering the new adjoint variables λi,i=1,2,…,6, into formula (22). They are not stored as the control variables ui, but as temporary vectors u¯i. The control variables ui are set as the convex combination of the last iteration of ui, namely, oldui, and the temporary vectors u¯i. That is,(23)ui=cku¯i+1−ckoldui,where k is the current iteration and c∈0,1. Here, we adopt c=0.2.The control periodT 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 A1=0.4,A2=0.1,B1=0.4,B2=0.4,B3=0.1,andB4=0.1. The numerical settings of Aii=1,2 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 Bii=1,2,3,4.We next consider three control schemes:(i)
Strategy 1. Only considering to reduce contact rates β1t and β2t, 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 β1t and β2t 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. The time-varying optimal control parameters u1, u2, u3, and u4 of Strategy 3 are shown in Figure 9.Figure 6
The time-series plots ofI (a) and L (b) without (red line) and with (blue line) control Strategy 1.
(a)
(b)Figure 7
The time-series plots ofI (a) and L (b) without (red line) and with (blue line) control Strategy 2.
(a)
(b)Figure 8
The time-series plots ofI (a) and L (b) without (red line) and with (blue line) control Strategy 3.
(a)
(b)Figure 9
The optimal control variablesui,i=1,2,3,4, of Strategy 3.Figure 10
The relation between the rate of vaccination withR0.


## 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 R0 and show that it can be used to govern the stochastic dynamics of the HFMD model (1): if R0<1, the system has a unique stable disease-free equilibrium which means the extinction of HFMD; if R0>1, 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 numberR0, 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. 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 [15–17].(iii)
The control strategy: we investigate the HFMD dynamics with three types of control schemes: (1) to reduce the contact ratesβ1t and β2t (see Figure 6); (2) to increase the recovery rates δ1 and δ2 (see Figure 7); (3) to reduce β1t and β2t and increase δ1 and δ2 (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 (β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 Figure3, 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.



---

*Source: 5902698-2020-02-24.xml*