Modelling and analysis of HFMD with the effects of vaccination, contaminated environments and quarantine in mainland China

1.   Introduction

Hand, foot, and mouth disease (HFMD) is a common childhood illness, mainly caused by coxsackievirus A16 (CVA16) and enterovirus 71 (EV71), as well as certain enteroviruses, including coxsackievirus A4, 5, 9, 10, B2 and 5 ^[1,2]. HFMD usually occurs in children under six years old, although it may also occur in older children and adults ^[3]. The usual incubation period is 2 to 7 days, and an infected individual will fully recover after 2 to 10 days. HFMD was first reported in New Zealand in 1957 and first named in America in 1959 ^[4], then it has widely spread in America ^[5], Europe ^[6,7] and Asia ^{[8,9,10,11,12,13,14,15,16,17]}.

Recently, HFMD has seriously threatened the health of infants and young children in mainland China. It has been classified as a class III infectious disease in the National Stationary Notifiable Communicable Diseases, and the annual data of HFMD have been archived by the Chinese Center for Disease Control and Prevention (CCDC) ^[3] and the National Health and Family Planning Commission of the People's Republic of China (NHFPC) ^[18] since 2008. Figure 1 shows the reported annual cumulative cases of more than 1.5 million since 2010 and more than 2 million since 2014. In analyses of analyze the epidemiological characteristics and pathogenic spectrum of FHMD in mainland China, EV71 and CVA16 have been considered major pathogens in the past decade, and the positive rates of EV71 and CVA16 were about 35%$\thicksim$55% and 20%$\thicksim$35%, respectively ^{[19,20,21,22,23]}. As described in ^[22,23], the EV71 vaccine was developed and has been used in mainland China since 2016, which have protected many young children protect away from the EV71 virus infections. And, from the data shown in Figure 1, we conclude that vaccination measures may have led to a lower number of cases in 2016 or 2017 than in 2014 or 2015, in contrast to the diseases trend in the past. The clinical trial results ^[22,23] showed that the vaccine provides immunity only for the EV71 virus; it provides no form of protection against infection for other viruses. Furthermore, there are few studies on the effects of EV71 vaccination on the spread of HFMD. It should be noted that epidemiological modelling is an important tool helpful to understand infectious disease spread and control ^[24,25,26].

A number of epidemiological models that do not consider EV71 vaccination have been established to investigate the transmission dynamics and to predict infections of HFMD. For instance, Wang and Sung (2006) ^[27] used a susceptible-infectious-recovered ($SIR$) model to analyse epidemic situations of FHMD in Taiwan from 1999 to 2003; Tiing and Labadin (2008) ^[28] studied $SIR$ model to predict the number of the infected and the duration of an outbreak in Sarawak Malaysia; Roy and Halder (2012) ^[29] proposed a susceptible-exposed-infectious-quarantine-recovered ($SEIQR$) model of HFMD; Ma and Liu (2013) ^[30] extended the $SEIQR$ model, considered asymptomatic infectious, to formulate a more realistic $SEII_{e}R$ model, which has been fitted to data of HFMD in Shandong Province, China; Wang and Xiao (2016) ^[31,32] further considered indirect transmission coming from the contaminated environment to establish an $SEII_{e}RW$ model. Moreover, some studies on the dynamical behaviors or application of the above models are shown in ^{[33,34,35,36,37,38,39,40]}. However, most of these models have considered either quarantine for infectious individuals ^[30,35,37] or contaminated environment caused by infectious individuals ^[31] alone. Moreover, no previous studies have considered the shared effects of vaccination, quarantine measures, and maintaining sanitation on disease. From the point view of practical application, the lack of key factors may affect the modeling.

To quantify these issues, we propose a new mathematical model, which extends the existing ones such as $SEIQR$ ^[30] and $SEII_{e}RW$ ^[31], to analyse the effects of EV71 vaccination, clean-up of contaminated environments, and quarantine on HFMD in mainland China simultaneously. In order to analyse the effects of EV71 vaccination, we take the individuals withEV71 vaccination as an independent compartment in modelling, which is a well deal with the case that EV71 vaccine against the EV71 virus but not against others. A combination of analytical and numerical techniques are used to analyze the proposed model. Finally, we use our model to individually fit the monthly reported data of HFMD in mainland China from 2015 to 2017. This paper is organized as follows: In Section 2, a novel epidemiological model of HFMD is proposed. In Section 3, the dynamic behaviors of this model are analyzed. In Section 4, the model is used to fit the real reported data, and some control strategies are discussed. In Section 5, we discuss and summarize our conclusions. In Appendix A, a new estimating algorithm based on BP neural network is designed.

2.   Model formulation

In this section, we propose an eight-compartment HFMD model, named $SVEII_{e}QRW$. The underlying structure of the model comprises classes of individuals are susceptible $S(t)$, vaccinated $V(t)$ (individuals with EV71 vaccine), exposed but not yet infectious $E(t)$, infectious with symptoms (symptomatic infectious) $I(t)$, infectious but not yet symptomatic (asymptomatic infectious) $I_{e}(t)$, infectious and hospitalized or quarantined or isolated $Q(t)$, and recovered $R(t)$. Using the same description in ^[31], $W(t)$ is the density of pathogen of the contaminated environment including door handles, towels, handkerchiefs, toys, utensils, bedding, underclothes, etc. at time $t$. Note that the individuals $V(t)$ vaccinated with EV71 vaccine can also be infected with FHMD by FHMD viruses except EV71.

A flow diagram describing the model is shown in Figure 2, and the model can be represented by the following ordinary differential equations:

where $k, p$ are nonnegative, other parameters are positive and its biological meanings are listed in Table 1.

Remark 1. Only the authors in ^[32] have considered the vaccination measure for HFMD, but it assumed that vaccinees against all HFMD viruses, in which the individuals being vaccinated in the susceptible class straight enter to the recovered class in modelling. However, it should be noted that the current vaccine is only resistant to the EV71 virus ^[22,23], which is well modeled in our model by adding a vaccinated class $V(t)$. In addition, the previous studies did not have taken into account the effects of both contaminated environment and quarantine on the spread of disease.

3.   Stability analysis and persistence

In this section, we discuss the dynamic behaviors for system (2.1). We first introduce some notations which will be used throughout this paper. $\mathbb{R}^{n}$ denotes the $n$-dimensional Euclidean space, and $\mathbb{R}^{n}_{+}\triangleq\{(x_{1}, x_{2}, \cdots, x_{n}) \in\mathbb{R}^{n}:x_{i}\geq 0, i = 1, \cdots, n\}$. $A = (a_{ij})_{n\times n}$ denotes a matrix of $n\times n$ real matrices. The superscript $T$ denotes matrix or vector transposition. For a continuous and bounded function $f:[0, +\infty)\rightarrow R$, denote $f_{\infty} = \liminf\limits_{t \to +\infty}f(t)$ and $f^{\infty} = \limsup\limits_{t \to +\infty}f(t)$. Let

be the initial condition of system (2.1). Clearly, any solution of system (2.1) with non-negative initial conditions is non-negative, and the following lemma shows that the solutions are uniformly ultimately bounded.

Lemma 3.1. The solutions of system (2.1) eventually enter

Proof. It is obvious that the population size $N(t) = S(t)+V(t)+E(t)+I(t)+I_{e}(t)+Q(t)+R(t)\geq 0$ and $W(t)\geq 0$. From system (2.1), the population size satisfies the following equation

According to the comparison theorem, there exists $t_{1}>0$ such that $N(t) \leq \frac{\Lambda}{d}$, for $t\geq t_{1}$.

It follows from the last equation of (2.1) that for $t\geq t_{1}$,

Similarly, there exists $t_{2}>t_{1}$ such that $W(t) \leq \frac{(\lambda_{1}+\lambda_{2})\Lambda}{d\delta}$, for $t\geq t_{2}$. Therefore, the solutions of system (2.1) are uniformly ultimately bounded. This completes the proof.

It is clear that system (2.1) has a disease-free equilibrium

Using the next generation matrix method developed by Van den and Watmough ^[41], we obtain the basic reproduction number

where

and

$R_{I}$ is the average number of secondary infected individuals generated by a symptomatic infected individual; $R_{I_{e}}$ is the average number of secondary infected individuals generated by an asymptomatic infected individual; $R_{3}$ is the average number of secondary infected individuals generated by free-living viruses in the environment shed by the infected individuals. $R_{0}$ is a central concept in measuring the transmission of infectious diseases.

In this following, we discuss the stability of the disease-free equilibrium and the persistence for system (2.1).

Theorem 3.2. The disease-free equilibrium $E_{0}$ is locally asymptotically stable if $R_{0} < 1$ and unstable if $R_{0}>1$.

Proof. The Jacobian matrix $J(E_{0})$ of system (2.1) at $E_{0}$ is given by

We can obtain the eigenvalues for (3.4) as $-(k+d)$, $-d$, $-(\gamma_{2}+d+m)$, $-(\eta_{1}+\eta_{2}+d)$, and roots of

where

Note that if $R_{0} < 1$, then $R_{1} < 1$, $R_{2} < 1$, $R_{3} < 1$. Hence, it can get that $a_{i}>0, i = 1, \cdots, 4$. Further, one has

According to Hurwitz criterion, all the roots of (3.5) have negative real part. Therefore, $E_{0}$ is locally asymptotically stable. If $R_{0}>1$, then $a_{4} < 0$. This means that (3.5) has at least one positive real part root. Then, $E_{0}$ is unstable. This completes the proof.

In order to analyse the global asymptotic stability of the disease-free equilibrium $E_{0}$ for system (2.1), denote

where $\hat{\beta}_{1} = \max(\beta_{1}, \bar{\beta}_{1}), \hat{\beta}_{2} = \max(\beta_{2}, \bar{\beta}_{2}), \hat{\sigma} = \max(\sigma, \bar{\sigma})$. It is obvious that $R_{0}\leq\hat{R}_{0}$. Moreover, it follows that $R_{0} = \hat{R}_{0}$ if and only if $\beta_{1} = \bar{\beta}_{1}, \beta_{2} = \bar{\beta}_{2}, \sigma = \bar{\sigma}$.

Theorem 3.3. If $\hat{R}_{0} < 1$, then the disease-free equilibrium $E_{0}$ is globally asymptotically stable.

Proof. From Theorem 3.2, we have that $E_{0}$ is locally asymptotically stable. Then, it only needs to prove the global attraction of the disease-free equilibrium $E_{0}$.

Note that any solution of system (2.1) with nonnegative initial is nonnegative and uniformly ultimataly bounded. By the fluctuation lemma ^[42], , there exists a sequence $\{t_{n}\}$ such that $t_{n}\rightarrow +\infty, S(t_{n})\rightarrow S^{\infty}$ and $\frac{dS(t_{n})}{dt}\rightarrow 0$ as $n\rightarrow +\infty$. It follows from the first equation of (2.1) that

Letting $n\rightarrow +\infty$, one has from (3.6) that

Similarly, it follows from the other equations of (2.1) that

By (3.7) and (3.8), one has

We claim that $E^{\infty} = 0$. Suppose not, it follows from (3.10), (3.11) and (3.14) that

Substituting inequalities (3.17), (3.18), (3.19) into (3.9), after some simplification, we derive that

By (3.20), one has that

It follows from (3.21) that

which is a contradiction since $S^{\infty}+V^{\infty}\leq\frac{\Lambda}{d}$ and $\hat{R}_{0} < 1$. This proves the claim. Since $E^{\infty} = 0$, then one has from (3.10)-(3.14) that $I^{\infty} = 0, I_{e}^{\infty} = 0, Q^{\infty} = 0, R^{\infty} = 0, W^{\infty} = 0$. Therefore,

By (3.15) and (3.16), one has

Moreover, using the fluctuation lemma again, there exists a sequence $\{s_{n}\}$ such that $s_{n}\rightarrow +\infty, S(s_{n})\rightarrow +\infty$, and $\frac{dS(s_{n})}{dt}\rightarrow 0$ as $n\rightarrow +\infty$. By the first and second equations of (2.1), one has

Letting $n\rightarrow +\infty$ and by using (3.22), we obtain that

It follows from (3.23) and (3.24) that

Thus,

This completes the proof.

Theorem 3.4. If $R_{0}>1$, then system (2.1) is uniformly persistent. That is, if $R_{0}>1$, there exists a small positive constant $\varepsilon>0$ such that $I_{\infty}>\varepsilon$, $I_{e\infty}>\varepsilon$ for system (2.1) with initial value $\varphi(0)$ and $I(0)>0$, $I_{e}(0)>0$.

Proof. In order to prove this result, the uniform persistence theorem in ^[44] is used. Define

Now we prove that system (2.1) is uniformly persistent with respect to $(X, X_{0})$.

First, it is easy to verify that both $X$ and $X_{0}$ are positively invariant for system (2.1), and $X_{0}$ is relatively closed in $X$. Moreover, by Theorem 3.1, the system (2.1) is point dissipative. Thus, the system (2.1) must exist a globally attractor. Set

Now we prove that

It is obvious that $M'_{\partial}\subseteq M_{\partial}$, then we only need to prove $M_{\partial}\subseteq M'_{\partial}$. Suppose not, let $\varphi(t)$ be a solution of system (2.1) with the initial condition $\varphi(0)$. Then, for any

and $\varphi(t)\neq M'_{\partial}$, there must exist at least one of $E(t)$, $I(t)$, $I_{e}(t)$, $Q(t)$, $R(t)$, $W(t)$ which is not zero. Without loss of generality, assume that $E(t) = 0$, $I(t) = 0$, $I_{e}(t) = 0$, $Q(t) = 0$, $R(t) = 0$, but $W(t)>0$. From system (2.1), one has the following equations, for $\forall t>0$,

This means that $\varphi(t) \not\in \partial X_{0}$ for $t>0$, which contradicts the assumption that $\varphi(t)\in M_{\partial}$. Hence, we obtain that $M_{\partial}\subseteq M'_{\partial}$. Since $M_{\partial} = M'_{\partial}$, we obtain that $M_{\partial}$ only has the disease-free equilibrium $E_{0}(S^{0}, V^{0}, 0, 0, 0, 0, 0, 0)$ and $E_{0}$ is compact and isolate invariant for $\varphi(0)\in M_{\partial}$.

Next, we only need to prove that $W^{s}(E_{0})\cap X_{0} = \emptyset$, where $W^{s}(E_{0})$ denotes the stable manifold of $E_{0}$. That is, there exists a positive constant $\varepsilon$ such that for any solution $\Phi_{t}(\varphi(0))$ of system (2.1) with the initial condition $\varphi(0)\in X_{0}$, one has

where $D$ is a distance function in $X_{0}$. Suppose not, then $D(\Phi_{t}(\varphi(0)), E_{0})^{\infty} < \bar{\varepsilon}$, for $\forall \bar{\varepsilon}>0$. Hence, there must exists $T>0$ such that $\frac{\Lambda}{k+d}-\bar{\varepsilon} \leq S(t)\leq \frac{\Lambda}{k+d}+\bar{\varepsilon}$, $\frac{k\Lambda}{d(k+d)}-\bar{\varepsilon} \leq V(t)\leq \frac{k\Lambda}{d(k+d)}+\bar{\varepsilon}$, $0\leq E(t)\leq \bar{\varepsilon}$, $0\leq I(t)\leq \bar{\varepsilon}$, $0\leq Q(t)\leq \bar{\varepsilon}$ and $0\leq W(t)\leq \bar{\varepsilon}$, for $t>T$. For $t>T$, we have

where

Consider an auxiliary system as follows

where $u = (u_{1}, u_{2}, u_{3}, u_{4})^{T}$ and

Recall that stability modulus of matrix $\bar{A}(\bar{\varepsilon})$, denoted by $S(\bar{A}(\bar{\varepsilon}))$, is defined by

Note that $\bar{A}(\bar{\varepsilon})$ is irreducible and has non-negative off-diagonal elements, then $S(\bar{A}(\bar{\varepsilon}))$ is a simple eigenvalue of $\bar{A}(\bar{\varepsilon})$ with a positive eigenvector (see Theorem A. 5 in ^[43]). By Lemma 2.1 in ^[44] and the proof of Theorem 2 in ^[41], the following equivalent inequations are hold:

Since $S(\bar{A}(\bar{\varepsilon}))$ is continous in $\bar{\varepsilon}$, choose $\bar{\varepsilon}>0$ small enough such that $S(\bar{A}(\bar{\varepsilon}))>0$ as $R_{0}>1$. Hence, let $u(t) = (u_{1}(t), u_{2}(t), u_{3}(t), u_{4}(t))^{T}$ be a positive solution of the auxiliary system (3.26), which is strictly increasing with $u_{i}(t)\rightarrow +\infty$ as $t\rightarrow +\infty$, $i = 1, \ldots, 4$. By the comparison principle, one has

By system (2.1), since

it is also gained that

This contracts with our assumption. Hence, $E_{0}$ is an isolated invariant set in $X$ and $W^{s}(E_{0})\cap X_{0} = \emptyset$.

Therefore, the system (2.1) is uniformly persistent if $R_{0}>1$. This completes the proof.

Remark 2. As far as the authors' knowledge, all the existing results of autonomous epidemiological models for HFMD such as ^{[28,29,31,34,35,40]} have obtained the similar result that the basic reproduction number $R_{0}$ is the only criterion to determine whether the disease is eliminated, i.e., the model only has a disease-free equilibrium if $R_{0} < 1$ and it is uniformly persistent if $R_{0}>1$. Because we consider vaccination and the vaccinated individuals can also be infected by HFMD viruses except EV71 in modeling, different from the existing results, the disease is extinct as $\hat{R}_{0} < 1$ (see Theorem 3.4). In addition, our model may have positive equilibriums even if the basic reproduction number $R_{0}$ is less than 1 (see the analyses in Section 5).

In particular, if we ignore the influence of vaccination, i.e., let $k = 0$ for system $(2.1)$ and define the model as $SEII_{e}QRW$, then the reproduction number for $SEII_{e}QRW$ is defined as

By using the similar analysis as that in the proof of 3.3 and Theorem 3.4, one can easily derive the following corollary. Its proof is omitted.

Corollary 1. If $\check{R}_{0} < 1$, then system $SEII_{e}QRW$ has an unique disease-free equilibrium $\check{E}_{0}(S^{0}, E^{0}, I^{0}, I^{0}_{e}, Q^{0}, R^{0}, W^{0}) = (\frac{\Lambda}{d}, 0, 0, 0, 0, 0, 0)$ and $\check{E}_{0}$ is globally asymptotically stable; If $\check{R}_{0}>1$, then system $SEII_{e}QRW$ is uniformly persistent.

From Corollary 1, it is clear that $\check{R}_{0}$ is the only threshold to decide whether the disease is extinct or not.

4.   Numerical analysis

In this section, by using system (2.1), we simulate the reported data of HFMD in mainland China. From the website of CCDC ^[3], and NHFPC ^[18], we have obtained the monthly numbers of newly reported HFMD cases in mainland China from January 2015 to December 2017. We only consider the population of young children under six years old, a high risk group for HFMD. Note that EV71 virus has been vaccinated in mainland China starting in 2016, then we assume that the vaccinated rate $k$ of system (2.1) is zero before 2016.

4.1.   Data fitting

First, based on BP neural network, we design an algorithm to estimate unknown parameters of system (2.1) (see Appendix A) to. Then, we estimate parameters of system (2.1) by minimizing the following error function

where $\hat{I}_{i}$ is the reported data of HFMD in the $i$ month. $I_{i}$ is numerically computed solutions of $I(t)$ of system (2.1) in the $i$ month. $N$ is the annual total population of less than 6 years old.

All the estimated parameters are shown in Table 2, and the error $ER$ for every year is computed and included in Table 3. It is obvious that $ER < 0.0001$ in every year. Therefore, our model fits well with the reported data of HFMD for each year from 2015 to 2017 in mainland China. Taking the year of 2017 as an example, the numerical simulation of the number of symptomatic infectious individuals $I(t)$ of HFMD by system (2.1) is shown in Figure 3.

4.2.   Sensitivity analysis

We all know that the basic reproduction number $R_{0}$ is an important quantity in characterizing the spread of disease.

Based on the basis of our parameter values, the basic reproduction number $R_{0}$ and basic components of $R_{0}$ include $R_{1}, R_{2}, R_{3}$ are calculated in every year in Table 4. Note that the basic reproduction $R_{0} = R_{1}+R_{2}+R_{3}$ for each year is greater than 1, where $R_{1}, R_{2}$, and $R_{3}$ with respect to symptomatic infected, asymptomatic infected individual, and free-living viruses in the environment, respectively. Moreover, in Table 4, the values of $R_{1}$ and $R_{2}$ are far greater than $R_{3}$, which means that the new patient has a high probability of being infected by symptomatic infectious individuals and asymptomatic infectious individuals, rather than the contaminated environment caused by infectious individuals.

Next, by using Latin hypercube sampling (LHS) and partial rank correlation coefficients (PRCCS) ^[46] (Marino et al., 2008), we investigate the effects of parameters on $R_{0}$. It notes that PRCCS, showing which parameters have the largest influence on model outcomes, is calculated using the rank transformed LHS matrix and output matrix ^[46]. Take 2000 simulations per round, choose a uniform distribution for all the parameters with ranges listed in Table 5, and test for significant PRCCs for all parameters of system (2.1).

Figure 4 shows that the PRCC values illustrate the dependence of $R_{0}$ on each input parameter. Assume that absolute values of PRCC$>0.4$ is highly negatively correlated between input parameters and $R_{0}$, then it is easy to know that $\beta_{1}, \beta_{2}, \sigma, \bar{\beta}_{1}, \bar{\beta}_{2}, \lambda_{2}$ are highly positively correlated to $R_{0}$, while $k, \gamma_{1}, \gamma_{3}, \delta, p$ are highly negatively correlated to $R_{0}$. That is, reducing positively correlated parameters and increasing negatively correlated parameters can reduce value of $R_{0}$, in which lower value of $R_{0}$ means lower spread of the epidemics.

4.3.   Prevention and control strategies

According to sensitivity analysis of $R_{0}$ in Subsection 4.2, we have known the key factors affecting the spread of HFMD. Then, some detailed numerical experiments provide effective measures for the prevention and control of the spread of FHMD in mainland China. Moreover, we pay more attention to that the varieties of vaccinated individuals, quarantined individuals, and environmental health impact HFMD epidemic in mainland China, which can be shown by adjusting the related parameters $k, p, \delta$ for system (2.1). Taking the year of 2017 as an example, the prevention and control strategies are proposed as follows:

Strategy 1: Increasing the vaccinated rate $k$ for the susceptible individuals $S(t)$ can reduce the spread of HFMD. Though the vaccinated individuals $V(t)$ can also be infected, but they have a lower infected probability than the unvaccinated individuals $S(t)$. In Table 2, all transmission rates of $V(t)$ of system (2.1) such as $\bar{\beta}_{1}, \bar{\beta}_{2}$ and $\bar{\sigma}$ are smaller than the transmission rates of $S(t)$. Figure 5(a) shows that the simulation of $I(t)$ of system (2.1) decreases with increasing the vaccinated rate $k$, while the solution curve of symptomatic infected individuals $I(t)$ of system (2.1) is at the top in the figure as $k = 0$. Moreover, Figure 5(b) and Figure 5(c) show that the peak values and annual cumulative cases of symptomatic infected individuals $I(t)$ decreases with the increase of $k$. Except for parameter $k$, Assuming that all the parameter of system (2.1) are shown in Table 1, and comparing Figure 5(d) with Figure 5(a), Figure 5(b) and Figure 5(c), we obtain the following results: (ⅰ) The size of $R_{0}$ determines the degree of the disease outbreaks, which are the same as the analyses in Subsection 4.2; (ⅱ) If the vaccinated rate $k < 0.07$, i.e., $R_{0} < 1$, then the disease persists (see Theorem 3.3); (ⅲ) If the vaccinated rate $k>0.07$, which means that $R_{0} < 1$ and $\hat{R}_{0}>1$, then we can not be sure whether the disease disappears or persists (see Theorem 3.3 and Theorem 3.4, and the reasons are analyzed in detail in Section 5). Furthermore, In order to investigate the effects of current EV71 vaccination strategy for HFMD in mainland China, we compare the annual simulated cumulative cases by system (2.1) without vaccination with the report cases. In Figure 6, the annual reported cumulative cases are lower than the annual simulation cumulative cases, and we also obtain that the current vaccine measures have reduced the total number of patients in 2016 and 2017 by 17% and 22%, respectively.

Strategy 2: Increasing the quarantined rate $p$ for the symptomatic infectious individuals $I(t)$ can reduce the spread of HFMD. From system (2.1), the quarantined individuals $Q(t)$ can not spread diseases and the symptomatic infectious individuals $I(t)$ are highly infectious, then transforming the symptomatic infectious into the quarantined is an efficient control strategy. From Figure 7(d), we obtain that the disease persists if the quarantined rate $p < 0.73$, while the long-term behaviors of the disease may not be ensure if $p>0.73$. The other analyses for Figure 7 are the same as strategy 1.

Strategy 3: Increasing the clearance rate $\delta$ of the pathogens caused by the infectious individuals containing $I(t)$ and $I_{e}(t)$ can reduce the spread of HFMD (see Figure 8). For instance, we can improve health-care education such as washing hands after using the toilet and before meals, and making air fresh indoors and so on. It should popularize health knowledge and advocate good personal hygiene habits in kindergardens, schools, hospitals and other places. Kindergardens should clean and disinfect toys and appliances every day. Moreover, hospitals should strengthen infection control practices to avoid nosocomial cross infection. In addition, from Figure 8(d), the disease persists.

In short, if we use the above prevention and control measures, the HFMD would be effectively controlled and the number of infections would decrease rapidly in short time. Those measures can effectively prevent and control the large-scale diffusion of HFMD in mainland China.

5.   Discussion and conclusion

In this paper, we have proposed an eight-compartment dynamic model for HFMD in mainland China. The proposed model extends the existing models by considering vaccination, contaminated environment and quarantine simultaneously, where the vaccinated individuals with EV71 vaccine can also be infected by other FHMD viruses and the transmission rates with respect to the vaccinated individuals may be different from the unvaccinated individuals. Hence, unlike the previous results, we take the vaccinated individuals as an independent compartment in modeling. The main purpose of this study is to examine the effect of vaccination, contaminated environment and quarantine on the spread of HFMD in mainland China.

By investigating the dynamic behaviors of our model, we have identified two critical parameters as the basic reproduction number $R_{0}$ and the threshold parameter $\hat{R}_{0}$ ($R_{0}\leq \hat{R}_{0}$). Then, the dynamic behaviors of system (2.1) could be divided into three cases (see Section 3):

(ⅰ) If $R_{0}\leq\hat{R}_{0} < 1$, then system (2.1) has only one disease-free equilibrium $E_{0}$ and $E_{0}$ is globally asymptotically stable, which means that the disease is extinct (see Theorem 3.3);

(ⅱ) If $1 < R_{0}\leq\hat{R}_{0}$, then system (2.1) is uniformly persistent (see Theorem 3.4);

(ⅲ) If $R_{0} < 1 < \hat{R}_{0}$, then dynamic behaviors of system (2.1) are difficult to determine.

Unlike the previous results of epidemiological models for HFMD ^{[28,29,31,34,35,40]}, in which the disease-free equilibrium is globally asymptotically stable if the basic reproduction number $R_{0} < 1$, in case (iii) that $E_{0}$ of system (2.1) may not be globally asymptotically stable and system (2.1) may have positive equilibriums. For example, if take $\Lambda = 1.3\times 10^{4}, \beta_{1} = 2.29\times 10^{-7}, \beta_{2} = 1.5\times 10^{-7}, \sigma = 1\times 10^{-7}, \bar{\beta}_{1} = 0.1\times 10^{-7}, \bar{\beta}_{2} = 1\times 10^{-7}, \bar{\sigma} = 0.7\times 10^{-7}, \alpha = 0.4, \rho = 0.5, \eta_{1} = 2.29, \eta_{2} = 1.5, \gamma_{1} = 0.9, \gamma_{2} = 2, \gamma_{3} = 0.9, \lambda_{1} = 0.04, \lambda_{2} = 0.04, \delta = 0.4, d = 1\times 10^{-3}, p = 0.1, k = 0.1, m = 1\times 10^{-4}$ for system (2.1), then we obtain a disease-free equilibrium

and a positive equilibrium (see Figure 9)

Based on the above parameter values, we calculate $R_{0} = 0.8981$ and $\hat{R}_{0} = 2.6993$. Then, in this case, when $R_{0} < 1 < \hat{R}_{0}$, it is obvious that system (2.1) has a disease-free equilibrium and a positive equilibrium, and none of them is globally asymptotically stable. The reason our model has dynamic behaviors different from those seen in the existing research is that our model considers vaccination, and the vaccinated individuals can also be infected by HFMD viruses except EV71; thus, the reproduction number R0 can not be the only threshold to judge whether the disease is extinct or not.

The proposed model has been used to individually fit the annual reported cases of HFMD in mainland China from 2015 to 2017. Moreover, the basic reproduction number $R_{0}$ of each year is greater than 1 (see Table 4), which means system (1) is uniformly persistent. Using Latin hypercube sampling (LHS) and partial rank correlation coefficients (PRCCS) for sensitive analysis of $R_{0}$, we find the key factors affecting transmission of HFMD, where the highly positively correlated parameters are $\beta_{1}, \beta_{2}, \sigma, \bar{\beta}_{1}, \bar{\beta}_{2}, \lambda_{2}$ and the highly negatively correlated parameters are $k, \gamma_{1}, \gamma_{3}, \delta, p$ (see Subsection 4.2). We analyse effect of the varieties of the rates $k, \delta$ and $p$, with respect to vaccinated individuals, quarantined individuals and environmental health respectively, on the spread of HFMD in mainland China. As a result, we deduce that increasing the rates $k, \delta, p$ can reduce the number of infectious individuals in detail in Subsection 4.3. Moreover, it also easily verifies that for any $k>0$ or $\delta>0$ or $p>0$, the threshold parameter $\hat{R}_{0}>1$. This means that our control strategies can reduce the spread of FHMD in mainland China, but the disease may not go extinct in mainland China. To eliminate the disease, it is necessary to reduce the rates $\beta_{1}, \beta_{2}, \sigma, \bar{\beta}_{1}, \bar{\beta}_{2}$ and $\gamma_{1}, \gamma_{3}$, but this is usually very difficult, in that the average contact rate of children may be hard to control and because in mainland China, medical conditions are difficult to improve over a short period.

It should be mentioned that we fitted our proposed model to the reported data and estimated some unknown parameters in each year, so that it is not necessary to consider seasonal infection variations of this disease. Note that HFMD is a typical seasonal disease, which is usually affected by temperature and humidity ^{[19,30,33,36]}. In the future, we will model periodic infection of the disease together with vaccination, contaminated environment, and quarantine. The authors thank the anonymous referees, whose careful reading, insights, valuable comments, and suggestions significantly enabled us to improve the quality of the paper.

Acknowledgments

H. Zhao is supported by the National Natural Science Foundation of China (No 11571170). D. Wu is supported by the Natural Science Foundation of Anhui Province of China (No 1608085MA14), and the Programs of Educational Commission of Anhui Province of China (No KJ2015A152). The authors thank anonymous reviewers for their valuable comments. Their comments have enhanced the clarity and the quality.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Appendix A.    

The BP neural network algorithm for estimating parameters

The parameter estimation methods for nonlinear problem that are most often used in practical applications are the nonlinear least squares estimator, Metropolis-Hastings algorithm, and so on ^[47]. The nonlinear least squares estimator is a parameter estimation method of nonlinear model parameters based on the criterion of error squared and minimum ^[48]. The Metropolis-Hastings algorithm estimates the posterior distributions of the parameters and the latent variables by using Markov chain Monte Carlo methods, which usually assumes that the sample satisfies certain probability distribution such as normal distribution, Gaussian distribution, chi-square distribution and so on ^[49]. Possibly due to limited by sample size, prior distribution assumptions or other constraints, the above two algorithms are used to estimate the unknown parameters of system (2.1) are not well, such as between the reported HFMD cases in mainland China and the simulation of $I(t)$ of system (2.1) in 2017, the fitting errors $ER$ with respect to the nonlinear least squares estimator (use the Fmincon function in Matlab) and metropolis-Hastings algorithm (use the source code in ^[47]) are $4.03\times10^{-4}$ and $7.46\times10^{-5}$, respectively, which are much larger than that of our estimated algorithm based on BP neural network (see Table 3). Because the neural network is not limited by the hypothesis of sample size and prior probability distribution, we design an algorithm based on BP neural network ^[50], which is one of the most widely neural network used and highly efficient ^[50], to estimate unknown parameters of system (2.1).

The topology of the BP neural network is shown in Figure 10. The BP neural network consists of three layers: input layer, hidden layer and output layer. The number of nodes in the input layer, the hidden layer and the output layer are 12, 93 and 8, respectively. For simplicity, we assume that $\theta = (\theta_{1}, \theta_{2}, ..., \theta_{m})^{T}$ denotes the vector of unknown parameters of system (2.1), where each element value must be set in reasonable interval. For example, set $\theta = (\beta_{1}, \beta_{2}, \sigma, \bar{\beta}_{1}, \bar{\beta}_{2}, \bar{\sigma}, \rho, k)$ in 2017 seen in Table 2. Let $I = (I_{1}, I_{2}, ..., I_{12})^{T}$ and $\hat{I} = (\hat{I}_{1}, \hat{I}_{2}, ..., \hat{I}_{12})^{T}$. Taking the year of 2017 as an example, the estimation algorithm is described as follows:

${\bf Step 1}$ Generate sample data: Generate $n$ group sample data $\theta^{i} = (\theta_{1}^{i}, \theta_{2}^{i}, ..., \theta_{m}^{i})^{T}$, $i = 1, 2, ..., n$ by Latin hypercube sampling. Then, taking $\theta^{i}$ into system (2.1), we can also obtain $n$ group numerically computed data $I^{i} = (I_{1}^{i}, I_{2}^{i}, ..., I_{12}^{i})^{T}$, $i = 1, 2, ..., n$ seen in Table 2;

${\bf Step 2}$ Use BP neural network to train data: Set $\bar{P} = (I^{1}, I^{2}, ..., I^{n})$ and $\bar{T} = (\theta^{1}, \theta^{2}, ..., \theta^{n})$, and each element of $\bar{P}$ and $\bar{T}$ is normalized into interval $[-1, 1]$. Then, taking $\bar{P}$ as input vector and $\bar{T}$ as output vector, a three-layers BP neural network is used to train dataset, and let $N_{et}$ be the training completed network, which the training error is less than the setting value. The detailed process of BP neural network can be seen in ^[50], and its specific settings are described below;

${\bf Step 3}$ : Estimate parameters: First, taking reported data $\hat{I} = (\hat{I}_{1}, \hat{I}_{2}, ..., \hat{I}_{12})^{T}$ of HFMD as the input vector into the training completed network $N_{et}$, we can gain a output forecast vector $\theta^{*} = (\theta_{1}^{*}, \theta_{2}^{*}, ..., \theta_{m}^{*})^{T}$. Then, taking $\theta^{*}$ into system (2.1), it get the simulation of $I_{i}$. Finally, we set a small positive content $\varepsilon = 0.0001$ to estimate the error function $ER$ (4.1). if $ER < \varepsilon$, the estimated parameters $\theta^{*}$ are reasonable, otherwise it goes to step 1 and reset the settings of BP neural network.

Figure 11 shows $n$ group numerically computed symptomatic data $I^{i}$ of system (2.1), $i = 1, 2, ..., n$, which are generated in step 1. In step 2, by using Matlab 2014b toolbox, the main settings of BP neural network are described as: (1) The transfer functions of hidden layer and output layer are tansig and purelin, respectively; (2) The network training function is trainlm; (3) The weight/bias learning function is learngdm; (4) The training error function is mse (mean squared error), and its value is 0.01; (5) The maximum iterations is 5000. Figure 12 shows the iteration process of the BP neural network, which the network complete training data before the setting maximum iterations 5000.