# Modeling Periodic HFMD with the Effect of Vaccination in Mainland China

Authors: Lei Shi; Hongyong Zhao; Daiyong Wu
Journal: Complexity
 (2020)

---

## Abstract

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


---

## Body

## 1. Introduction

Hand, foot, and mouth disease (HFMD) is a common infectious disease among infants and children caused by intestinal viruses of the Picornaviridae family. There exist many types of HFMD virus, such as coxsackievirus A5, A10, A16, A19 types and EV71 type, and so on. Among these viruses, the main viruses causing HFMD are coxsackievirus A16 (CVA16) and enterovirus 71 (EV71) [1]. HFMD mostly occurs in nursery schools or kindergartens, and its high incidence seasons are summer and autumn [1]. The usual incubation period is 2 to 7 days, and the recovery period is 7 to 10 days [2]. Since HFMD was first reported in New Zealand in 1957 [3], it has spread around the world, especially in Asia [4–7].HFMD is a mild, self-limiting illness that primarily affects infants and young children. In recent years, HFMD is also a common illness in mainland China. Since 2008, the government has classified it as a class III infectious disease in the National Stationary Notifiable Communicable Diseases (NSNCD) [8], and the monthly cases have been archived by the Chinese Center for Disease Control and Prevention (CCDC) [9] (see Figure 1). The main viruses causing HFMD in mainland China are CVA16 and EV71, and the proportions of HFMD cases caused by CVA16 and EV71 are about 30% and 44%, respectively [10–12]. EV71 vaccine, which efficacy against EV71-associated HFMD was 97.4%, has been developed and put into use in mainland China beginning in 2016, so many young children would be protected from EV71 infections [13, 14]. As of May 2019, monovalent EV71 vaccines have not yet been included in the routine paediatric vaccination programme in China, which means that they are needed to be paid out-of-pocket by parents [15–17]. Meanwhile, vaccine coverage of these monovalent EV71 vaccines among children aged 6 months to 5 years ranges from <10% to 50% in different provinces in mainland China [15–17].Figure 1
The confirmed monthly cumulative HFMD cases from the Chinese Center for Disease Control and Prevention (CCDC) from January 2008 to June 2019.Compartmental differential equation modeling is an important tool for understanding infectious disease spread and control. There are two types of dynamic models for HFMD. One is an autonomous ordinary differential model where all the parameters such as transmission rate, birth population, recovery rate, and so on of whom are assumed as constants. For instance, Tiing and Labadin [2] used a simple SIR model to predict the number of the infected and the duration of an outbreak in Sarawak Malaysia; Roy and Halder [18, 19] proposed a deterministic susceptible-exposed-infectious-recovered (SEIR) model and a susceptible-exposed-infectious-quarantine-recovered (SEIQR) model of HFMD, and Viriyapong and Wichaino [20] further analyzed dynamic behaviors of SEIR; Yang et al. [21] and Li et al. [22] researched the optimal control of SEIQR and its dynamic behavior, respectively; Li et al. [23] further constructed a two-stage-structured SEIQR model to fit the HFMD data in mainland China from 2009 to 2014; Wang et al. [24] further considered indirect transmission coming from the contaminated environments to establish a SEIIeRW model (Ie is the number of the asymptomatic infectious individuals), and Sharma and Samanta [25] considered the combined effect of asymptomatic infectious individuals and quarantine measures on the spread of HFMD. It should be noted that the reported monthly cases of HFMD from CCDC [9] exhibit the periodic pattern. From a practical point of view, clarifying the mechanisms that link seasonal environmental changes to disease dynamics may aid in forecasting the long-term health risks, crucial in developing an effective public health program and in setting objectives, and utilizing limited resources more effectively. For these reasons, the other approach is a nonautonomous differential equation model, which considers periodic transmission rates. Of these, Liu [26] constructed a periodic SEIQR model to simulate the dynamics of HFMD transmission and showed that quarantine has a positive impact on controlling the spread of HFMD; Zhu et al. [27] used the periodic SEIQR model to investigate the spread of seasonal HFMD in Wenzhou city, Zhejiang Province, China; Ma et al. [28], considering asymptomatic infectious individuals (Ie), proposed a periodic SEIIeQR model to analyse HFMD in Shandong Province, China, which shows that asymptomatic infectious subpopulation plays an important role in the spread of HFMD. However, among the two types of dynamic models, fewer studies considered the influence of EV71 vaccination on the spread of HFMD.Recently, Samanta [29] and Wang et al. [30] considered vaccination in modeling for HFMD. In these studies, the vaccinated individuals were assumed to be immunized against all HFMD-associated enteroviruses, in which a vaccinated individual would be transferred from the susceptible compartment to the recovered compartment in modeling (see system (1) in [29] or system (1) in [30]). However, in mainland China, only EV71 vaccine existed, and the vaccine only has well immunity for EV71, but not for others [13, 14, 31, 32]. Thus, while there are multiple HFMD-associated enteroviruses, only one vaccine for EV71 is available. A natural problem is how to build a suitable model to reflect the influence of EV71 vaccination on the spread of HFMD with multiple pathogenic viruses in mainland China.From the practical point of view, periodic transmission functions are needed to consider in modeling. From Figure1, it should be noted that the monthly reported cases of HFMD in mainland China from CCDC [9] exhibit periodic patterns. The annual cycle patterns or seasonal patterns may be attributed to these facts such as climatic, geographical, and demographic information but not limited to them [28, 30, 33]. Although different seasonality across the whole country, Xing et al.’s epidemiological study of HFMD in China, 2008–2012 [33], implied that geographical differences in seasonal patterns were weakly associated with climate and demographic factors (variance explained 8–23% and 3–19%, respectively). Our analysis of previous studies for periodic HFMD models [4, 26, 28, 30] finds that all the chosen periodic transmission functions have a similar form, and its period is 1 year, which effectively can reflect the annual seasonal changes for HFMD. However, annual cycle patterns aside, the epidemics of HFMD in mainland China have displayed the following two phenomena: (1) the data from CCDC [8] (see Figure 1) show that HFMD has a two-year epidemic cycle because the epidemic period of partial HFMD-associated enteroviruses in some areas is two years [34, 35]; (2) Figure 1 shows that the number of cases suddenly rises in September and October each year because of the opening and closing of schools in mainland China [36]. The previous models based on the disease’s transmission functions [4, 26, 28, 30] cannot simulate the above two phenomena. Hence, to better reveal the mechanism of the epidemic, it is very important to find a more suitable periodic transmission function for modeling.Given the issues discussed above, we propose a new mathematical model to investigate the effect of EV71 vaccination on the spread of HFMD in mainland China. To analyze the effect of EV71 vaccination, we derive our modeling ideas from the previous modeling for influenza [37]. Then, in modeling, we divide the total individuals into two groups, the unvaccinated and the vaccinated. In addition, novel two-year periodic transmission functions, with considering the effects of school opening and closing, are also considered in our modeling. Then, we evaluate the basic reproduction number R0 to analyze the dynamical behaviors of the model and use the proposed model to fit the monthly reported data of HFMD in mainland China from January 2008 to June 2019. Furthermore, through sensitivity analysis of the basic reproduction number R0 in terms of key parameters, we explore some effective prevention and control measures for HFMD in mainland China.The paper is organized as follows: in Section2, we introduce a new periodic epidemiological model of HFMD. In Section 3, the dynamic behaviors of this model are analyzed theoretically. In Section 4, simulations of the model, sensitivity analysis of the basic reproduction number, and some prevention and control measures are performed. In Section 5, we discuss and summarize our conclusions.


## 2. Model Formulation

Recently, Samanta [29], Wang et al. [30] considered the effect of the vaccination in modeling for HFMD, in which it assumed that vaccine can resist all the HFMD-associated enteroviruses. However, there only exists EV71 vaccine to be used in mainland China, and the vaccine only has well immunity for EV71 but not for others [13, 14]. It notes that the modeling techniques in [29, 30] cannot be used to really reflect the spread of the disease with a single effective vaccine and multiple viruses. In order to model the influence of EV71 vaccination on the spread of HFMD with multiple pathogenic viruses in mainland China, we divide the total individuals into two different groups, the unvaccinated and the vaccinated, in modeling. Next, based on the above analyses, we denote the total number of unvaccinated and vaccinated individuals by Nst and Nvt, respectively, and classify each of them into four subclasses: susceptible, exposed, infectious with symptoms (symptomatic infectious), and infectious but not yet symptomatic (asymptomatic infectious), with the number of unvaccinated individuals denoted by St,Est,Ist,Iset and Rst and the number of vaccinated individuals denoted by Vt,Evt,Ivt,Ivet and Rvt, respectively. Denote the total population size Nt=Nst+Nvt.In addition, many epidemiological models for HFMD ([26–28, 30]), considering annual seasonal changes, were simulated by using the sinusoidal function β1t=a01+a1sinπ/6t with one year period (take a month as the unit time). As discussion in Section 1, seasonal changes and school holidays will change the contact rate of children, and the virus contagion changes periodically, some of which are one year for a cycle and some of which are two years for a cycle. Thus, taking a month as the unit time, we define periodic transmission functions with the period of 2 years as β1t=a01+a1sinπ/2t+a2sinπ/6t+a3sinπ/12t+a4sin−1+π/3t and β2t=b01+b1sin−1+π/2+π/4t+b2sinπ/6t+b3sinπ/12t for symptomatic infectious individuals and asymptomatic infectious individuals, respectively. The studies in [31] showed that, in mainland China, the EV71 vaccine was targeted at infants or children aged 6 to 71 months. They were also a high-risk group for HFMD [31, 32]. Next, let Λ be the monthly number of children entering the age of 6–71 months, and d be progression rate leaving the children group aged 6–71 months. We also assume that p is the vaccinated rate for infants or children aged 6–71 months. Moreover, the results in [15, 17] showed that vaccine coverage of these monovalent EV71 vaccines among children aged 6 months to 5 years ranges from <10% to 50% in different provinces in mainland China in recent years. This implies that the current vaccination rate of the monovalent EV71 vaccine is very low in mainland China. In the case of the low vaccine coverage, we let ks be the proportion of HFMD cases in mainland China, without vaccination, caused by HFMD-associated enteroviruses except EV71. Furthermore, we assume that a symptomatic or an asymptomatic infectious individual with vaccination has the same recovery rate as a symptomatic or asymptomatic infectious individual without vaccination. Then, the transmission dynamics associated with these subpopulations are shown in Figure 2.Figure 2
Flow diagram representing transmission routes and other processes modeled by system (1).The model is described by the following ordinary differential equations:(1)dSdt=1−pΛ−β1tSIs+Iv−β2tSIse+Ive−dS+ηsRs,dEsdt=β1tSIs+Iv+β2tSIse+Ive−d+αEs,dIsdt=αρEs−d+μ+λ1Is,dIsedt=α1−ρEs−d+λ2Ise,dRsdt=λ1Is+λ2Ise−d+ηsRs,dVdt=pΛ−β1tVksIs+Iv−β2tVksIse+Ive−dV+ηvRv,dEvdt=β1tVksIs+Iv+β2tVksIse+Ive−d+αEv,dIvdt=αρEv−d+μ+λ1Iv,dIvedt=α1−ρEv−d+λ2Ive,dRvdt=λ1Iv+λ2Ive−d+ηvRv,where p is nonnegative and other parameters are positive, and its biological meanings are listed in Table 1.Table 1
Descriptions and values of parameters in system (1).
Para.
Value
Interpretation (units)
Source
Λ
1.41
×
10
6
The monthly number of children entering the age of 6–71 months (/month)
[38]
D
1.30
×
10
−
2
Progression rate leaving the children group aged 6–71 months (/month)
[38]
k
s
0.56
Proportion of infectious individuals, without vaccination, caused by HFMD-associated enteroviruses except EV71 (none)
[39]
μ
3.4
×
10
−
4
Disease-related death rate of HFMD (/month)
[8, 9]
p
0.23
Vaccinated rate for children aged 6–71 months (/month)
Fitting
η
1
3.50
×
10
−
2
Rate fromRst to St (/month)
[24]
η
2
3.50
×
10
−
2
Rate fromRvt to Vt (/month)
Assumption
1
/
α
4/30
Average incubation period (/month)
[2]
ρ
0.23
Proportion of symptomatic infectious individuals (none)
Fitting
λ
1
3.53
Recovery rate of symptomatic infectious individuals containingIst and Ivt (/month)
[2]
λ
2
3.46
Recovery rate of asymptomatic infectious individuals containingIset and Ivet (/month)
Fitting
a
0
1.09
×
10
−
7
The baseline contact rate of symptomatic infectious individuals containingIst and Ivt (none)
Fitting
a
1
0.12
The magnitude of forcing of symptomatic infectious individuals containingIst and Ivt (none)
Fitting
a
2
0.22
The magnitude of forcing of symptomatic infectious individuals containingIst and Ivt (none)
Fitting
a
3
0.03
The magnitude of forcing of symptomatic infectious individuals containingIst and Ivt (none)
Fitting
a
4
0.15
The magnitude of forcing of symptomatic infectious individuals containingIst and Ivt (none)
Fitting
b
0
2.01
×
10
−
8
The baseline contact rate of asymptomatic infectious individuals containingIset and Ivet (none)
Fitting
b
1
0.08
The magnitude of forcing of asymptomatic infectious individuals containingIset and Ivet (none)
Fitting
b
2
0.19
The magnitude of forcing of asymptomatic infectious individuals containingIset and Ivet (none)
Fitting
b
3
0.13
The magnitude of forcing of asymptomatic infectious individuals containingIset and Ivet (none)
Fitting


## 3. Stability Analysis and Persistence

In this section, we will investigate the dynamic behaviors of system (1). First, we introduce some notations which will be used throughout this paper. ℝn denotes the n-dimensional Euclidean space, and ℝ+n≜x1,x2,…,xn∈ℝn:xi≥0,i=1,…,n. A=aijn×n denotes a n×n real matrix. The superscript T denotes matrix or vector transposition. Let(2)φ0=S0,Es0,Is0,Ise0,Rs0,V0,Ev0,Iv0,Ive0,Rv0T∈ℝ+10,be the initial condition of system (1), and(3)φt=St,Est,Ist,Iset,Rst,Vt,Evt,Ivt,Ivet,RvtT∈ℝ+10,be the solution of system (1) at time t for t≥0.From the first equation of system (1), if St=0, then dS/dt≥0. According to function continuity and monotonicity, when S0≥0, we have St≥0. Moreover, we have similar results for other state variables of system (1). Therefore, any solution of system (1) with nonnegative initial conditions is nonnegative. Next, the following lemma shows that the solutions are uniformly ultimately bounded.Lemma 1.
The solutions of system (1) eventually enter(4)Γ=St,Est,Ist,Iset,Rst,Vt,Evt,Ivt,Ivet,RvtT∈ℝ+10:0≤St+Est+Ist+Iset+Rst≤1−pΛd,0≤Vt+Evt+Ivt+Ivet+Rvt≤pΛd.
Moreover,Γ is a positively invariant set.Proof.
It is obvious that the unvaccinated population sizeNst=St+Est+Ist+Iset+Rst≥0 and the vaccinated population size Nvt=Vt+Evt+Ivt+Ivet+Rvt≥0. From system (1), we have(5)dNsdt=1−pΛ−dS−dEs−d+μIs−dIse−dRs=1−pΛ−dNs−μIs≤1−pΛ−dNs,(6)dNvdt=pΛ−dV−dEv−d+μIv−dIve−dRv=pΛ−dNv−μIv≤pΛ−dNv.
According to the comparison theorem, it follows from (5) and (6) that(7)Nst≤Ns0−1−pΛde−dt+Λd,(8)Nvt≤Nv0−pΛde−dt+Λd.
Using (7) and (8), we obtain that limsupt⟶∞Nst=1−pΛ/d and limsupt⟶∞Nvt=pΛ/d. Hence, the solutions of system (1) are uniformly ultimately bounded. Moreover, if the initial condition φ0∈Γ, then the solution φt∈Γ, i.e., Γ is a positively invariant set for system (1). This completes the proof.
It is obvious that system (1) always has a disease-free equilibrium:(9)E0S0,Es0,Is0,Ise0,Rs0,V0,Ev0,Iv0,Ive0,Rv0=1−pΛd,0,0,0,0,pΛd,0,0,0,0T.
Next, we first introduce a very important threshold basic reproduction numberR0. The basic reproduction number R0 is the number of secondary cases in which one case would produce in a completely susceptible population [40]. We calculate the basic reproduction number R0 for system (1) by using the general calculation procedure given by Wang and Zhao [41]. It is easy to verify that system (1) satisfies the conditions (A1)–(A7) given in [41], and we have(10)ℱ=β1tSIs+Iv+β2tSIse+Ive00β1tVksIs+Iv+β2tVksIse+Ive00,V=d+αEs−αρEs+d+μ+λ1Is−α1−ρEs+d+λ2Ised+αEv−αρEv+d+μ+λ1Iv−α1−ρEv+d+λ2Ive.
Thus, we obtain that(11)F¯=0β1tS0β2tS00β1tS0β2tS00000000000000β1tksV0β2tksV00β1tV0β2tV0000000000000,V¯=d+α00000−αρd+μ+λ10000−α1−ρ0d+λ2000000d+α00000−αρd+μ+λ10000−α1−ρ0d+λ2.
LetΦV¯t be the monodromy matrix of the linear ω-periodic system dz/dt=−V¯tz, and ρΦV¯ω be the spectral radius of ΦV¯ω, where ω is the period. Let Yt,s be the 6×6 matrix solution of the following initial value problem:(12)dYt,sdt=−V¯Yt,s,∀t≥s,Ys,s=I,where I is the 6×6 identity matrix.
Suppose thatϕs is the initial distribution of infectious individuals, then Fsϕs is the rate of new infections produced by the infectious individuals introduced at time s, and Yt,sFsϕs represents the distribution of those infectious individuals who are newly infected at time s and still remain in the infected compartment at time t for t≥s. Then, one has that(13)ψt≔∫−∞tYt,sFsϕsds=∫0∞Yt,t−aFt−aϕt−ada,is the distribution of accumulative new infections at time t produced by all those infected individuals ϕs introduced at time previous to t. Let Cω be the ordered Banach space of all ω periodic functions from ℝ6 to ℝ with maximum norm ⋅ and positive cone Cω+≔ϕt∈Cω:ϕt≥0,∀t∈ℝ. Then, we can introduce a linear operator L:Cω⟶Cω by(14)Lψt≔∫0∞Yt,t−aFt−aϕt−ada,∀t∈ℝ,ϕ∈Cω.
Following the results obtained in [41], the basic reproduction number R0 for model (1) is defined as the spectral radius of operator L, i.e., R0≔ρL.
In order to analyze the dynamic behaviors of system (1), we have the following result.Lemma 2 (Theorem 2.2 in [41]).
The following statements are valid:(i)
R0=1 if and only if ρΦF¯−V¯ω=1(ii)
R0>1 if and only if ρΦF¯−V¯ω>1(iii)
R0<1 if and only if ρΦF¯−V¯ω<1
Therefore, the disease-free equilibriumE0 of system (1) is locally asymptotically stable if R0<1 and unstable if R0>1. Next, we will analyze the global dynamics of the disease-free equilibrium E0.Theorem 1.
IfR0<1, then the disease-free equilibrium E0 of system (1) is globally asymptotically stable.Proof.
According to Lemma2, one has that R0<1 is equivalent to ρΦF¯−V¯ω<1. Hence, we only need to prove that when ρΦF¯−V¯ω<1, E0 is globally asymptotically stable.
Chooseε>0 small enough such that ρΦF¯−V¯+εMωω<1, where(15)Mω=0β1tβ2t0β1tβ2t0000000000000β1tksβ2tks0β1tβ2t000000000000.
From Lemma1, for ∀ε>0, there exists t1>0 such that for t>t1, St≤1−pΛ/d+ε and Vt≤pΛ/d+ε. Thus, for t>t1, it follows from system (1) that(16)dEsdt≤β1t1−pΛd+εIs+Iv+β2t1−pΛd+εIse+Ive−d+αEs,dIsdt=αρEs−d+μ+λ1Is,dIsedt=α1−ρEs−d+λ2Ise,dEvdt≤β1tpΛd+εksIs+Iv+β2tpΛd+εksIse+Ive−d+αEv,dIvdt=αρEv−d+μ+λ1Iv,dIvedt=α1−ρEv−d+λ2Ive.
Consider the following auxiliary system:(17)dudt=Aεh,where h=h1,…,h6T,Aε=F¯−V¯+εMω. Following Lemma 2.1 in [42], there exists a positive ω-periodic function h^t=h^1,…,h^6T such that ht=eμth^t is a solution of system (17), where μ=1/ωlnρΦF¯−V¯+εMωω. Due to the continuity of functions and following Lemma 2, when ε⟶0 and R0<1, then ρΦF¯−V¯+εMωω<1. Hence, we have ut⟶0 as t⟶∞, which implies that the zero solution of system (17) is globally asymptotically stable. Based on the comparison principle, it follows from system (17) that(18)limt⟶∞Est,Ist,Iset,Evt,Ivt,IvetT=0.
Using (18), it follows from the fifth and last equations in system (1) that(19)limt⟶∞Rst=0,limt⟶∞Rvt=0.
Using (18) and (19), it follows from the first and sixth equations in system (1) that(20)limt⟶∞St=1−pΛd,limt⟶∞Vt=pΛd.
Thus,E0 is globally asymptotically stable when R0<1. This completes the proof.
Next, we will analyze the persistence of system (1). Define(21)X0=S,Es,Is,Ise,Rs,V,Ev,Iv,Ive,RvT∈Γ:Es>0,Is>0,Ise>0,Rs>0,Ev>0,Iv>0,Ive>0,Rv>0,∂X0=Γ\X0.
Assume thatut,φ0 is the solution of system (1) with the initial value φ0. According to the fundamental existence-uniqueness theorem in [43], ut,φ0 is unique. Denote that P:Γ⟶Γ is the Poincaré map with respect to system (1), i.e.,(22)Pφ0=uω,φ0,∀φ0∈Γ,where ω is the period. It is obvious that(23)PmE0=umω,E0,∀m≥0.
It follows from Lemma1 that Γ is positively invariant and P is point dissipative for system (1). In order to analyze the persistence for system (1), we give the following lemma.Lemma 3.
IfR0>1 , then there exists δ>0 such that for any φ0∈X0 with φ0−E0≤δ, we have(24)limsupm⟶∞DPmφ0,E0≥δ,where D is a distance function in X0.Proof.
According to Lemma1, if R0>1, then ρΦF¯−V¯ω>1. Choose ϵ>0 small enough such that ρΦF¯−V¯−εMωω>1. Now, we prove that(25)limsupm⟶∞DPmφ0,E0≥δ.
If it is false, then(26)limsupm⟶∞DPmφ0,E0<δ,for some φ0∈X0. Assume that there exists m>0 such that DPmφ0,E0<δ. According to the continuity of the solutions associated with the initial conditions, when φ0−E0≤δ, it follows that(27)ut1,Pmφ0−ut1,E0<ε,m≥0,t1∈0,ω.
For∀t≥0, let t=mω+t1, where t1∈0,ω and m=t/ω, which is the greatest integer less than or equal to t/ω. Then, for ∀t≥0, we have(28)ut,φ0−ut,E0=ut1,Pmφ0−ut1,E0<ϵ,t1∈0,ω,which implies that 1−pΛ/d−ϵ≤St≤1−pΛ/d+ϵ,pΛ/d−ϵ≤Vt≤pΛ/d+ϵ,t≥0. Then, for φ0−E0≤δ, one has(29)dEsdt≥β1t1−pΛd−εIs+Iv+β2t1−pΛd−εIse+Ive−d+αEs,dIsdt=αρEs−d+μ+λ1Is,dIsedt=α1−ρEs−d+λ2Ise,dEvdt≥β1tpΛd−εksIs+Iv+β2tpΛd−εksIse+Ive−d+αEv,dIvdt=αρEv−d+μ+λ1Iv,dIvedt=α1−ρEv−d+λ2Ive.
Consider the following auxiliary linear system:(30)dgdt=A¯ϵg.
As the similar proof process in Theorem 1, there exists a positive ω-periodic functiong^t=g^1,…,g^6T such that gt=eμtg^t is a solution of system (30), where μ=1/ωlnρΦF¯−V¯−εMωω. Since μ=1/ωlnρΦF¯−V¯−εMωω>1, it follows that if g0>0, then gt⟶∞ as t⟶∞. From the comparison principle, we have(31)Est⟶∞,Ist⟶∞,Iset⟶∞,Evt⟶∞,Ivt⟶∞,Ivet⟶∞,as(32)Es0>0,Is0>0,Ise0>0,Ev0>0,Iv0>0,Ive0>0,for t⟶∞, which is a contradiction with Lemma 1. Therefore, limsupm⟶∞DPmφ0,E0≥δ. This completes the proof.Theorem 2.
IfR0>1, then system (1) is uniformly persistent and has at least one positive periodic solution.Proof.
First, we prove that system (1) is uniformly persistent, i.e., there exists δ>0 such that any solution φt of system (1) with the initial value φ0∈X0 satisfies(33)lim inft⟶∞Est,Ist,Iset,Rst,Evt,Ivt,Ivet,RvtT≥δ,δ,δ,δ,δ,δ,δ,δT.
For∀φ0∈X0, solving system (1), we obtain that(34)St=e−∫0td+β1tIsu+Ivu+β2tIseu+IveuduS0+∫0t1−pΛ+ηsRse∫0td+β1tIsu+Ivu+β2tIseu+Iveududu,Est=e−α+dtEs0+∫0tβ1tSuIsu+Ivu+β2tSuIseu+Iveudu,Ist=e−d+μ+λ1tIs0+∫0tαρEsudu,Iset=e−d+λ2tIse0+∫0tα1−ρEsudu,Rst=e−d+ηstRs0+∫0tλ1Isu+λ2Iseudu,Vt=e−∫0td+β1tksIsu+Ivu+β2tksIseu+IveuduV0+∫0tpΛ+ηsRse∫0td+β1tksIsu+Ivu+β2tksIseu+Iveududu,Evt=e−α+dtEv0+∫0tβ1tVuksIsu+Ivu+β2tVuksIseu+Iveudu,Ivt=e−d+μ+λ1tIv0+∫0tαρEvudu,Ivet=e−d+λ2tIve0+∫0tα1−ρEvudu,Rvt=e−d+ηstRv0+∫0tλ1Ivu+λ2Iveudu.
From (34), it is obvious that for t≥0, the solution φt∈X0 for system (1) as the initial value φ0∈X0. This means that X0 is positively invariant. In addition, ∂X0 is relatively closed in Γ, and following Lemma 1, the discrete-time system Pmm≥0 admits a global attractor in Γ. Then, we set(35)M∂≔φ0∈∂X0:Pmφ0∈∂X0,∀m≥0.
Now, we prove that(36)M∂=S,0,0,0,0,V,0,0,0,0T∈∂X:S≥0,V≥0≜M∂′.
It is easy to know thatM∂′⊆M∂, then we only need to prove M∂⊆M∂′. That is, for any(37)φ0=S0,Es0,Is0,Ise0,Rs0,V0,Ev0,Iv0,Ive0,Rv0T∈∂X0,we have(38)Ismω=Isemω=Rsmω=Evmω=Ivmω=Ivemω=Rvmω=0.
If not, then there existsm1≥0 such that(39)Ism1ω,Isem1ω,Rsm1ω,Evm1ω,Ivm1ω,Ivem1ω,Rvm1ωT≥0.
If we replace the initial timet=0 with t=m1ω, then it follows from (34) that(40)Est>0,Ist>0,Iset>0,Rst>0,Evt>0,Ivt>0,Ivet>0,Rvt>0,for ∀t>m1ω. This means that Pmφ0∉∂X0, for ∀m>m1, which contradicts the definition of M∂. Therefore, it follows that M∂⊆M∂′. Moreover, there only exists one fixed point E0 of P in M∂.
According to Lemma3, E0 is an isolated invariant set in Γ. From the acyclicity theorem on uniform persistence for maps (Theorem 1.3.1 in [44]), it follows that P is uniformly persistent with respect to X0,∂X0. Moreover, according to Theorem 1.3.1 in [44], the solutions of system (1) are uniformly persistent with respect to X0,∂X0, i.e., there exists δ>0 such that any solution φt of system (1) with the initial value φ0∈X0 satisfies(41)liminft⟶∞Est,Ist,Iset,Rst,Evt,Ivt,Ivet,RvtT≥δ,δ,δ,δ,δ,δ,δ,δT.
Next, we prove that system (1) has a positive periodic solution. Theorem 1.3.6 in [44] implies that P has a fixed point(42)S∗0,Es∗0,Is∗0,Ise∗0,Rs∗0,V∗0,Ev∗0,Iv∗0,Ive∗0,Rv∗0T∈X0.
Moreover, we can also obtain thatS∗0>0 for system (1). Suppose not, assume S∗0=0. By the periodicity of P, we have S∗0=S∗nω,n=1,2,3,…. From the first equation of (34), one has(43)S∗nω=e−∫0nωd+β1tIsu+Ivu+β2tIseu+IveuduS∗0+∫​0t1−pΛ+ηsRse∫0nωd+β1tIsu+Ivu+β2tIseu+Iveududu>0,which contradicts with S∗0=S∗nω=0. Thus, it follows that S∗0>0. Then, from (34), it is easy to obtain that S∗t>0,Es∗t>0,Is∗t>0,Ise∗t>0,Rs∗t>0,V∗t>0,Ev∗t>0,Iv∗t>0,Ive∗t>0,Rv∗t>0, for t≥0. Due to the definition of semiflow P, we obtain that(44)S∗t,Es∗t,Is∗t,Ise∗t,Rs∗t,V∗t,Ev∗t,Iv∗t,Ive∗t,Rv∗tT,is a positive ω-periodic solution of system (1). This completes the proof.


## 4. Simulations and Sensitivity Analysis

In this section, by using system (1), we simulate the reported data of HFMD in mainland China from January 2008 to June 2019, predict the trend of the disease, and seek for some prevention and control strategies. We obtain the monthly number of newly reported HFMD cases in mainland China from the website of the Chinese Center for Disease Control and Prevention (CCDC) [9] and the National Health and Family Planning Commission of the People’s Republic of China (NHFPC) [8]. We only consider the population of young children aged 6–71 months, who are a high-risk group for HFMD [31, 32] and get demographic information from the National Bureau of Statistics of China (NBSC) [38]. From NBSC [38], we obtain the monthly number of children entering the age of 6–71 months Λ=1.41×106 and the average monthly rate of leaving the children group aged 6–71 months d=1.30×10−2. From CCDC [9] and NHFPC [8], we get the disease-related death rate μ=3.40×10−4. The proportion of infectious individuals without vaccination caused by HFMD viruses except EV71 is ks=0.56 [39]. Other parameters can be obtained from the literature or assumed on the basis of common sense (see Table 1). Unfortunately, there is no public data on the vaccinated rate p and the proportion of symptomatic infectious individuals ρ. Therefore, we estimate the values of unknown parameters containing p,ρ and β1,β2. By using the BP neural network algorithm [45] (see Appendix), we estimate parameters of system (1) by calculating the following error function:(45)ER=1n∑i=1nIi−I^i2,where I^i is the reported data of HFMD at time t. Ii is the numerically computed solution of It of system (1) at time t. n is the total number of sample points collected at all times.Next, we set the initial conditions for system (1). We take January 2008 as the initial time. It notes that the EV71 vaccine has been vaccinated in mainland China starting in 2016, then vaccinated rate p of system (1) is zero before 2016 [13, 14]. Hence, since the total number of vaccinated individuals Nvt=0 before January 2016, then it follows that V0=0,Ev0=0,Iv0=0,Ive0=0,Rv0=0. The initial number of the susceptible human population aged 6–71 months of 2008 is calculated as S0=9.23×107 from NBSC [8], and the number of the initial unvaccinated symptomatic infectious humans Is0=2.10×104 is from CCDC [9]. However, the initial number of the unvaccinated exposed humans Ev0, the unvaccinated asymptomatic infectious humans Ise0, and the unvaccinated recovered humans Rs0 cannot be obtained. We derive Ev0 reversely by the parameter αρ; get Ise0 by using the proportion parameter ρ to compare with Is0, and Rs0 is estimated roughly. The numerical simulation of the model on the number of HFMD cases in mainland China from January 2008 to June 2019 is shown in Figure 3, and its fitting error ER (45) is calculated as 6.71×10−4. In addition, in order to test how well system (1) actually reflects the data, by using the method in [22, 46], we consider the following hypotheses.Figure 3
The comparison between the reported HFMD cases in mainland China and the simulation ofIst+Ivt of system (1) from January 2008 to June 2019. The values of parameters are given in Table 1. The initial conditions are S0=9.23×107,Es0=8×104,Is0=2×104,Ise0=6×104,Rs=1×105,V0=0,Ev0=0,Iv0=0,Ive0=0,Rv0=0. It should be noted that we take k=0 before 2016.Null hypothesis,H0: the estimated parameters of system (1) are equal to actual values.Alternative hypothesis,H1: the estimated parameters of system (1) are not equal to actual values.Since the EV71 vaccine was only started in 2016, we consider hypothesis test for the two time periods from January 2008 to December 2015 and January 2016 to June 2019, and the chi-square values and degrees of freedom are shown in Table2. Then, from Table 2, we obtain that it cannot reject the null hypothesis at the 5% significant level by Pearson’s criterion of chi-square test [47]. Therefore, it indicates that, with these parameter values, there is a good fit between the simulation of model (1) and the HFMD cases in mainland China. Furthermore, based on our parameters listed in Table 1, we calculate the basic reproduction number as R0=ΦF¯−V¯24=1.09 after 2015, for this case where the vaccination is considered, which means that system (1) admits at least a positive periodic solution (see Theorem 2).Table 2
Chi-square values and degrees of freedom from January 2008 to December 2015 and January 2016 to June 2019.
Year
Jan. 2008–Dec. 2015
Jan. 2016–June 2019
Chi-square value
68.07
33.19
DF
82
28
AV
104.14
41.34
Note: AV denotes the accepting value is at 5% significant level with degrees of freedom DF.It notes the initial conditions exceptIs0 adopted in model fitting are mostly assumed and back-extrapolated by parameters. Thus, it is necessary to analyze the influence of the initial values S0,Es0,Is0,Ise0,Rs0 on the number of infected cases Ist+Ivt, which is illustrated in Figure 4. Figure 4 shows that the initial value S0 has a greater impact on Ist+Ivt while other initial values have little or no impact on Ist+Ivt. This implies that the selection of the initial values S0 is very important for our simulation. According to [6, 7, 10, 11, 32], due to most of the infected individuals of HFMD are 6–71 months old, then we carefully select that the susceptible individuals are 6–71 months old.Figure 4
The influence of initial conditions of humans on the number of HFMD casesIst+Ivt: (a) different values of S0; (b) different values of Es0; (c) different values of Is0; (d) different values of Ise0; (e) different values of Rs0.
(a)
(b)
(c)
(d)
(e)In order to investigate the effect of the EV71 vaccination strategy on HFMD in mainland China, we compare the trends of HFMD between with vaccination and without vaccination by simulating model (1), which is presented in Figure 5. In Figure 5, the solution curve without considering vaccination is on top of the one with considering vaccination after 2015, which means that the EV71 vaccination strategy is effective in reducing the spread of HFMD in mainland China in the last two years. In addition, we calculate the basic reproduction number R0 of system (1) without considering vaccination (p=0) as 1.24 and with considering vaccination (p=0.23) as 1.09, so system (1) in each case admits a positive periodic solution according to Theorem 2 (see the solution curves from 2030 to 2035 in Figure 5). In addition, Figure 6 shows that under the current vaccination measures, the total number of patients in 2016, 2017, and 2018 has reduced about 4.62×105, 5.16×105, and 7.39×105, respectively.Figure 5
The tendency of HFMD from system (1) with vaccination (p=0 before 2016 and p=0.23 in others) and without vaccination in a long time (p=0 in all times). The values of other parameters and the initial conditions are the same as in Figure 3.Figure 6
The comparison between the annual simulation cumulative cases by system (1) without vaccination (p=0) and the annual reported cumulative cases for HFMD in mainland China from 2016 to 2018. The values of other parameters and the initial conditions are the same as in Figure 3.To compare the performance of the periodic transmission functions in this paper with some previous research studies, by using system (1), we consider the following two cases:Periodic functions (I): choose the transmission functions in this paper such asβ1t=a01+a1sinπ/2t+a2sinπ/6t+a3sinπ/12t+a4sin−1+π/3t and β2t=b01+b1sin−1+π/2+π/4t+b2sinπ/6t+b3sinπ/12t, where the values of ai,bi,i=1,…,3 are given in Table 1.Periodic functions (II): choose the transmission functions in the previous research studies [26–28, 30] such as β¯1t=a¯01+a¯1sinπ/6t and β¯2t=b¯01+b¯1sinπ/6t. Based on our parameter values in Table 1, we estimate that the optimal values of a¯0,a¯1,b¯0,b¯1 are 1.19×10−7,0.23,0.21,2.03×10−8, respectively.The studies in [48] showed that the trend of monthly new cases of infectious diseases is usually closely related to the transmission functions. It notes that the periods of β¯1t and β¯2t are 12 months, and these functions have only one peak in each year. Meanwhile, the periods of β1t and β2t are 24 months, and these functions have two peaks in each year. As explained in Introduction, the transmission functions β1t,β2t are closer to the trend of actual data change in mainland China than β¯1t,β¯2t. Moreover, in Figure 7, the simulation of Ist+Ivt of system (1) is based on the periodic functions (I) (i.e., β1t,β2t) but the periodic functions (II) (i.e., β¯1t,β¯2t) can well fit the two phenomena, where HFMD has a two-year epidemic cycle and the number of the cases per year has a slight upward trend in September and October in mainland China (see the reasons in Introduction).Figure 7
The tendency of HFMD from system (1)with periodic functions (I) and with periodic functions (II). The values of other parameters and the initial conditions are the same as in Figure 3.Moreover, the period of the solution for system (1) with periodic functions (II) is 12 months, which fails to reflect the 24 months periodic variation of the disease, while these biological phenomena are well presented by the solution for system (1) with periodic functions (I). Furthermore, we calculate the fitting error ER for system (1) with periodic functions (I) and with periodic functions (II) which is 4×10−3 and 1.20×10−2, respectively. On the whole, our model, especially based on our selection of the transmission rate functions, can well reveal the objective laws of the spread of HFMD in mainland China.According to Theorems1 and 2, we have known that the basic reproduction number R0, which determines the outbreak of the disease or extinction, is an important quantity in characterizing the spread of disease. Therefore, to find some effective prevention and control measures for HFMD in mainland China, we perform sensitivity analysis to determine the influence of some parameters on R0. Due to some parameters of system (1) for HFMD are difficult or not to control such as Λ,d,α,ρ and so on, biologically, then we only consider the influence of the parameters p,β1,β1,λ1 on R0, which are depicted in Figure 8, and its related prevention and control measures are shown in Figure 9. Under the assumption of low vaccine coverage, choosing p≤0.3, Figure 8(a) shows that R0 decreases with the increase of p. This implies that the larger vaccinated rate p (p≤0.3) is, the less number of symptomatic infectious individuals Is+Iv is (see Figure 9(a)). Moreover, Figure 8(a) also shows that R0 is always larger than 1 as p≤0.3, which means that without other prevention and control, the disease is not extinct. Figures 8(b) and 8(c) show that the less the baseline contact rate of β1 or β2 is, the less R0 is. Moreover, when the baseline contact rate a0<9×10−8 or b0<1.45×10−8, R0 can be less than 1, which means that the disease is extinct. However, based on our estimating parameter values, a0 or b0 is far greater than 9×10−8 or 1.45×10−8, which means that HFMD becomes an endemic disease in mainland China. If the relevant departments can effectively control the contact rate between the healthy children and the infectious children and improve hygienic precautions and environmental cleaning, such as washing hands before meals and after using the toilet, and making air fresh indoors and so on, the disease can be effectively control or extinct (see Figures 9(b) and 9(c)). Figure 8(d) shows that R0 decreases with the increasing of the recovery rate of symptomatic infectious individuals λ1, i.e., if average recovery period 1/λ1 of symptomatic infectious individuals is less, then R0 is less, and the disease is extinct if λ1>4.28. Figure 9(d) shows that the number of symptomatic infectious cases decreases as the recovery rate λ1 increases, especially the disease is extinct when λ1=4.73. Therefore, we suggest that the patients should seek medical advice in time to achieve diagnosis and treatment early. At the same time, because the disease does not have targeted drugs, some targeted drugs need to be developed actively. Finally, in order to further discuss the credibility of our results, we will evaluate and discuss the consistency of the scope of R0 with other studies in China. Although the ideas of dynamic modeling for different studies based on the HFMD cases in China are not the same, the range of R0 of models may have similar characteristics. For instance, the study in [30], which assumed that the vaccine was against all HFMD-associated enteroviruses and analyzed the HFMD cases in mainland China from 2010 to 2014, showed that as the vaccination rate, the baseline contact rate of symptomatic infectious individuals and the baseline contact rate of asymptomatic infectious individuals change, the ranges of R0 are 1.15,1.58, 1.37,1.58, and 0.48,1.58, respectively. The study in [28], based on the HFMD cases in Shandong Province, China, which mainly discussed the impacts of quarantine measures and the fraction of asymptomatic infectious individuals on the spread of HFMD, indicated that, as the quarantine rate and the fraction of asymptomatic infectious individuals change, the ranges of R0 are 1.04,1.14 and 1.02,1.05. In a word, most studies for HFMD in China (see paper [27, 28, 30, 49]) showed that the range of R0 is around 1 and not far away from 1, which is similar to the change of R0 shown in Figure 8. Therefore, these discussions can also indirectly improve the credibility of our results.Figure 8
The influence of parameters onR0 (a) versus p; (b) versus a0, being the baseline contact rate of β1; (c) versus b0, being the baseline contact rate of β2; (d) versus λ1. Other parameter values are given in Table 1.
(a)
(b)
(c)
(d)Figure 9
The influence of parameters on symptomatic infectious individualsIs+Iv (a) versus p; (b) versus a0, being the baseline contact rate of β1; (c) versus b0, being the baseline contact rate of β2; (d) versus λ1. Other parameter values are given in Table 1, and the initial conditions are the same as in Figure 3.
(a)
(b)
(c)
(d)


## 5. Discussion

In this paper, we have proposed and analyzed a HFMD model with periodic transmission rates to take into account seasonal outbreak of HFMD infection and to consider EV71 vaccination. The studies in [29, 30] have assumed that the vaccine is effective against all enteroviruses, which may not be in line with the current situation of HFMD transmission in mainland China. In fact, one vaccine is available in mainland China and it works only for EV71 but not for others. To investigate the effects of EV71 vaccination on HFMD in the presence of multiple pathogenic viruses, our proposed model divided the total individuals into the unvaccinated and the vaccinated, which were described in different compartments (see Section 2). Moreover, we have chosen novel periodic transmission functions, unlike the previous studies in [4, 26, 28, 30], to reflect the two-year periodic nature and the infections’ change with the opening and closing of schools in mainland China (see Sections 1 and 2).We have analyzed the global dynamics of system (1) by evaluating the basic reproduction number R0≔ρL: if R0<1, then the disease-free equilibrium of system (1) is globally asymptotically stable (see Theorem 1), which means that the disease is extinct; if R0>1, then system (1) uniformly persists and has at least one positive periodic solution (see Theorem 2), which means that the disease persists or has not been eliminated.By fitting our model to the reported data on symptomatic cases of HFMD in mainland China, we estimated the basic reproduction numberR0=1.24 from 2008 to 2015 and R0=1.09 from 2016 to June 2019. The major factors for this difference are that EV71 vaccination has been used since 2016, the reason of which can also be seen in the analyses for the influence of the vaccinated rate p on R0 in Figure 8(a). Moreover, the periodic transmission functions chosen in our paper perform efficiently for system (1) to fit the reported data in mainland China, given the analyzes of Figure 7. Furthermore, to find effective prevention and control measures for HFMD in mainland China, we did sensitivity analysis to determine the influence of some artificial control parameters on R0. It showed that under the assumption of low vaccine coverage, increasing the vaccinated rate with respect to EV71 vaccination can reduce the number of infections, but the disease does not go extinct. If we can effectively control the contact rate between the healthy children and the infectious children (see the analyses of Figures 8(b) and 8(c) and Figures 9(b) and 9(c)) and improve the medical conditions (see the analyses of Figures 8(d) and 9(d)), the disease may be eliminated. In a word, the proposed prevention and control measures may be useful for public health governance.It notes that the observed aggregated monthly HFMD cases in Figure1 are a mixture of different seasonality in mainland China. Also, Xing et al.’s study of HFMD in mainland China from 2008 to 2012 [33] showed that seasonal diversities (although this was not the only factor) made the spread of HFMD present different characteristics in different regions, such as the annual amplitude of HFMD epidemics increased with increasing latitude and semiannual periodicity was the strongest in the south, and so on. It is difficult for us to use the data in Figure 1 to discuss the differences in the spread of HFMD in different regions of the country. Thus, we used the national average transmission rate functions β1t,β2t in our modeling, where β1t,β2t are periodic functions. In future studies, in order to avoid the influence of seasonal differences in different regions on the research results, we will consider a specific region such as a province in China.



---

*Source: 8763126-2020-01-09.xml*