Effects of behaviour change on HFMD transmission

ABSTRACT

We propose a hand, foot and mouth disease (HFMD) transmission model for children with behaviour change and imperfect quarantine. The symptomatic and quarantined states obey constant behaviour change while others follow variable behaviour change depending on the numbers of new and recent infections. The basic reproduction number ${\mathcal{R}}_0$ of the model is defined and shown to be a threshold for disease persistence and eradication. Namely, the disease-free equilibrium is globally asymptotically stable if ${\mathcal{R}}_0\le 1$ whereas the disease persists and there is a unique endemic equilibrium otherwise. By fitting the model to weekly HFMD data of Shanghai in 2019, the reproduction number is estimated at 2.41. Sensitivity analysis for ${\mathcal{R}}_0$ shows that avoiding contagious contacts and implementing strict quarantine are essential to lower HFMD persistence. Numerical simulations suggest that strong behaviour change not only reduces the peak size and endemic level dramatically but also impairs the role of asymptomatic transmission.

KEYWORDS:

• HFMD
• behaviour change
• basic reproduction number
• infection size
• global stability

MATHEMATICS SUBJECT CLASSIFICATIONS:

• 37N25
• 34D20
• 92D30

1. Introduction

The emergence and spread of infectious diseases not only pose a great threat to public health, but also hinder economic growth and social stability. Among them is hand, foot and mouth disease (HFMD), a common viral infection in infants and children under 5 years old. More than 20 types of enteroviruses can cause HFMD, of which coxsackievirus A16 (CA16) and enterovirus 71 (EV71) are the two major etiological agents. The disease is primarily spread through close contact with infected people, through the respiratory tract or the faecal-oral route. The incubation period for HFMD is typically around 3–5 days. The virus can be isolated from the faeces or throat secretions of patients during the incubation period. Most infected people are asymptomatic or have mild symptoms, but they can still spread the infection. The clinical manifestations of HFMD usually include fever, mouth sores and skin rash on hands, feet and mouth, lasting 7-10 days. In particular, CA16 is associated with mild cases while EV71 infections could be severe and even fatal. Recovery from infection induces immunity against that specific serotype but is not immune to other serotypes. No specific treatment is yet available for HFMD and pain relievers can be given to reduce discomfort. So far three inactivated monovalent EV71 vaccines have been approved in China, but the coverage is relatively low. Good personal hygiene such as washing hands properly can reduce the risk of infection substantially.

HFMD was first reported in New Zealand in 1957. The disease is now endemic in many parts of the world. In particular, it is prevalent in Southeast Asian countries and regions such as Singapore, Vietnam, Mongolia and Malaysia. For example, an outbreak of HFMD in the state of Sarawak in Malaysia in 1997 caused 2626 cases and 31 deaths [38]. There are also cases of HFMD in developed economies including Japan, the United States and Canada [32]. In China, the earliest recorded HFMD case was from Shanghai in 1981, and the first HFMD outbreak occurred in Tianjin in 1983, which resulted in more than 7,000 confirmed cases of CA16 infection. Since then, cases have emerged in Northeast, East and South of China. In Spring 2008, a large-scale outbreak of HFMD due to EV71 occurred in Fuyang City, Anhui Province, with almost 25,000 reported cases. Subsequently, HFMD was categorized as a legally notifiable disease in China and there are usually more than one million reported cases every year [18] (see Figure 1).

Figure 1. Annual number of HFMD cases reported in mainland China, 2008–2021.

Mathematical modelling of HFMD has a relatively short history compared to other common infectious diseases like influenza, measles and malaria, but it has developed rapidly in the past two decades. In 2004, Wang and Sung [41] adopted an SIRS model with seasonally-forced transmission rate to evaluate the impact of weather on severe enterovirus infections in Taiwan from 2000 to 2003. Tiing and Labadin [38] used an SIRS model that consists of natural birth and death and disease-caused death to predict the number of infectives and the duration of the outbreak in Sarawak, Malaysia. Roy and Halder [33] established a HFMD model to distinguish between symptomatic cases and asymptomatic cases. Liu [23] proposed a periodic SEIQRS model and showed that quarantine is beneficial to the control of HFMD. Ma et al. [27] generalized the model of Liu [23] to an SEIAQRS model by adding asymptomatic infection and transmission and fitted the model to the weekly HFMD data of Shandong Province. Li et al. [22] constructed a two-stage-structured (i.e. children and adults) SEIQRS model and estimated the annual basic reproduction number of HFMD in China from 2009 to 2014. To explore the impact of contaminated environments on HFMD dynamics, Wang et al. [43] formulated an SEIARS-W model with both human-to-human and environment-to-human transmission routes where W denotes the concentration of pathogen in the environment. Shi et al. [34] took vaccination, contaminated environments, quarantine and asymptomatic infection into consideration and evaluated the contribution of EV71 vaccine in reducing HFMD cases. Recently, Zhao et al. [48] developed a stage-structured (scattered and school children) SEIRS model where transmission rates and latent periods are temperature-dependent or time-periodic and found that the time-averaged system may underestimate the infection risk and size.

Human behaviour and disease transmission interact with each other strongly [28]. During the 1918 influenza pandemic, a range of nonpharmaceutical interventions such as isolating infected people, closing schools, restricting or banning mass gatherings, and promoting cleaning and disinfection, were tried to mitigate the disease spread across the U.S. and shown to reduce the overall mortality significantly in cities where early and effective interventions were implemented [4]. In response to the SARS epidemic in 2003, residents in Hong Kong quickly adopted preventive measures including mask use, frequent hand washing, and avoiding crowded places [20]. Travel restrictions, stay-at-home orders, personal protection equipment, and social distancing have been the key measures against the spread of SARS-CoV-2 [31]. The inclusion of behavioural changes in epidemiological models has grown fast in recent years [12]. We highlight some typical deterministic models in which behavioural change is incorporated into contact rates since a similar approach will be used in the current study. As early as 1973, London and Yorke [25] explored the mechanism of periodic disease outbreaks by modelling seasonal variation in contact rate due to the opening and closing of school. To consider the psychological effects, Capasso and Serio [7] introduced an SIR epidemic model with a nonlinear force of infection depending on the instantaneous number of infectives. Hsu and Hsieh [17] proposed a SARS model with quarantine and other intervention measures where the contact rate decreases with respect to the cumulative number of probable cases representing behaviour change by individuals. Liu et al. [24] considered another type of SARS model in which the transmission coefficient is an exponential decreasing function of the reported numbers of the exposed, infectious and hospitalized individuals. They showed that the presence of media/psychological impact could induce multiple outbreaks or even sustained periodic oscillations. To evaluate the effects of information and education campaigns on the HIV epidemic in Uganda, Joshi et al. [19] divided the susceptible population into three categories according to behaviour changes and added a compartment to reflect the amount of education campaigns. Cui et al. [10] analysed an SIS model where the contact rate is decreasing in the number of infected individuals and found that contact reduction driven by media coverage cannot change disease dynamics but can reduce the endemic level. Gao and Ruan [16] considered the SIS model with media coverage in a patchy environment and showed that increasing media coverage in any patch will reduce the infection size in each patch. Brauer [5] proposed a simple SIR epidemic model with different constant contact reductions for susceptible and infectious members. Misra et al. [30] established a model that divides the susceptible class into aware and unaware groups and includes the time delay in execution of awareness programmes. Xiao et al. [45] proposed an SIR model where the media effect is described by a piecewise smooth function that depends on the number of cases and its change rate. Based on an SEIR model, Collinson and Heffernan [9] compared the effects of different media functions on epidemic outcomes. Wang et al. [42] incorporated behaviour change into a cholera model via disease contact rates and host shedding rate. The role of behaviour change in containing disease spread has received considerable attention in modelling COVID-19 transmission [15].

Behaviour changes are common in the public health response to HFMD in China. For example, during the high transmission season, more frequent disinfection and sanitation practices and more careful health checks are performed in preschools and kindergartens. Children are encouraged to wash their hands frequently and avoid cross-contamination before eating and after playing outdoors. In the case of a serious epidemic, public health, education and communication departments work closely and run health campaigns intensively and release situation reports regularly. When a case is identified on campus, the patient is required to be treated at hospital or stayed at home and the entire class may be quarantined. The aim of this paper is to investigate the influence of behaviour change on the spread of HFMD. Our paper is structured as follows. In Section 2, we develop a mathematical model with status-dependent behaviour change to describe HFMD transmission. Section 3 establishes the threshold dynamics of the model system in terms of the basic reproduction number. In Section 4, data fitting and numerical simulations are conducted to further analyse the consequence of different control strategies. We conclude the paper in Section 5 with a discussion of our main findings and their implications.

2. Model formulation

The incidence rates of HFMD vary significantly among different age groups [18]. The vast majority of HFMD reported cases are among infants and children. Thus, we only consider the population consisting of children under 12 in this study. Since the case fatality rate of HFMD is very low, especially after the distribution of EV71 vaccines, we ignore the disease-induced mortality. The total population, $N(t)$, is divided into six mutually exclusive classes $S(t),E(t),I_s(t),I_a(t),Q(t)$ and $R(t)$ denoting the number of susceptible, exposed, symptomatically infected, mildly or asymptomatically infected, quarantined and recovered individuals at time t, respectively. So, $N(t)=S(t)+E(t)+I_s(t)+I_a(t)+Q(t)+R(t).$A susceptible individual can get infected through contact with symptomatic, asymptomatic or quarantined cases. After the incubation period, an exposed individual may or may not develop symptoms. Symptomatic cases will be quarantined after being diagnosed. Quarantine can reduce transmission but cannot perfectly prevent it. People who recovered from infections are assumed to acquire permanent immunity due to the lower incidence rate in elder children. The HFMD transmission is illustrated in Figure 2.

Figure 2. Flow chart for the HFMD model.

Let $U_S,U_E,U_{I_s},U_{I_a},U_Q$ and $U_R$ denote the contact rate of susceptible, exposed, symptomatic, asymptomatic, quarantined and recovered individuals, respectively. The normal contact rate under no behaviour change is denoted by c. In reality, an individual's response to an epidemic is affected by the person's health status [5]. People without symptoms or ‘healthy people’ (including susceptible, exposed, asymptomatic and recovered classes) are worried about contracting the disease, so they tend to take more protective measures as the epidemic is getting worse. Meanwhile, symptomatic cases have a high probability of being confirmed and their behaviours may no longer change with the progression of the epidemic. Based on these observations, we assume that $\begin{array}{rl}{U}_S& =U_E=U_{I_a}=U_R=f(I_s,Q)c,\\ U_{I_s}& ={\alpha }_{1}c,\\ U_Q& ={\alpha }_{2}c,\end{array}$where $f(I_s,Q)$ is a continuously differentiable, positive and decreasing function of $I_s$ and Q, and ${\alpha }_1\in (0,1)$ and ${\alpha }_2\in (0,1)$ are the relative contact rate of the symptomatically infected and quarantined populations, respectively. In other words, symptomatically infected, quarantined and the remaining individuals decrease their contact rate by a fraction of ${\alpha }_1,{\alpha }_2$ and $f(I_s,Q)$, respectively. The level of behaviour change is governed by the severity of an epidemic which is closely related to the number of reported cases and its change rate [45]. Since asymptomatic infections are hard to detect, only symptomatic cases are assumed to be reported. The cumulative number of reported cases that are still infectious is Q and the number of newly reported cases per unit time is $\kappa I_s$. Hence, we suppose that the behaviour change function f for healthy people depends on $I_s$ and Q. In particular, $f(0,0)=1$, i.e. no behaviour change occurs if there is no symptomatic or quarantined case (no reported case). Denote the transmission probability from an infectious individual to a susceptible individual per contact by p. For convenience, we assume the symptomatic and asymptomatic cases have the same infectivity. The theoretical results throughout this paper remain valid if otherwise. Thus, the usual transmission coefficient is $\beta =cp$. The total number of contacts over the entire population per unit time is $\begin{array}{rl}\mathcal{U}& =U_{S}S+U_{E}E+U_{I_s}{I}_s+U_{I_a}{I}_a+U_{Q}Q+U_{R}R\\ & =f(I_s,Q)c(S+E+I_a+R)+{\alpha }_{1}c{I}_s+{\alpha }_{2}cQ\\ & =f(I_s,Q)cN+({\alpha }_1-f(I_s,Q))c{I}_s+({\alpha }_2-f(I_s,Q))cQ.\end{array}$Since the sum of $I_s$ and Q only accounts for a small proportion of the total population, the approximation $U_S/\mathcal{U}\approx 1/N$ holds. Accordingly, the forces of infection driven by symptomatically infected, asymptomatically infected and quarantined classes are $\begin{array}{rl}{U}_S\frac{U_{I_s}{I}_s}{\mathcal{U}}p& \approx {\alpha }_{1}c\frac{I_s}{N}p=\frac{{\alpha }_1\beta I_s}{N},\\ U_S\frac{U_{I_a}{I}_a}{\mathcal{U}}p& \approx f(I_s,Q)c\frac{I_a}{N}p=\frac{f(I_s,Q)\beta I_a}{N},\\ U_S\frac{U_{Q}Q}{\mathcal{U}}p& \approx {\alpha }_{2}c\frac{Q}{N}p=\frac{{\alpha }_2\beta Q}{N},\end{array}$respectively. The above derivation and the flow chart lead to the following endemic model with nonnegative initial conditions for HFMD transmission (1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q,R)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q,R)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& ={\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q-\mu R,\end{array}$(1) where $\lambda (S,E,I_s,I_a,Q,R)=\beta \frac{{\alpha }_1{I}_s+f(I_s,Q)I_a+{\alpha }_{2}Q}{N}$is the force of infection for the susceptible population. All model parameters except $f(I_s,Q)$ are constants with their definitions and ranges summarized in Table 1.

Adding all equations of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q,R)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q,R)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& ={\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q-\mu R,\end{array}$(1) ) gives $\frac{\mathrm{d}N}{\mathrm{d}t}=\mathrm{\Lambda }-\mu N.$Since $N(t)=\frac{\mathrm{\Lambda }}{\mu }+(N(0)-\frac{\mathrm{\Lambda }}{\mu }){\mathrm{e}}^{-\mu t}\to N^{\ast }:=\frac{\mathrm{\Lambda }}{\mu }$ as $t\to \mathrm{\infty }$, it follows from the theory of asymptotically autonomous systems (see e.g. Castillo-Chavez and Thieme [8]) that system (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q,R)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q,R)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& ={\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q-\mu R,\end{array}$(1) ) can be reduced to (2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) where $\lambda (S,E,I_s,I_a,Q)=\beta \frac{{\alpha }_1{I}_s+f(I_s,Q)I_a+{\alpha }_{2}Q}{N^{\ast }}$ and $N^{\ast }=\frac{\mathrm{\Lambda }}{\mu }$. We can easily prove that system (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) is mathematically well-posed and biologically meaningful.

Theorem 2.1

For any initial condition starting in $\mathrm{\Omega }=\{(S,E,I_s,I_a,Q)\in {\mathbb{R}}_{+}^5|S+E+I_s+I_a+Q\le \mathrm{\Lambda }/\mu \}$, system (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) has a unique solution that remains in Ω for all time $t\ge 0$.

Table 1. Descriptions and ranges of model parameters (time unit is day).

	Description	Range	References
Λ	Recruitment rate	300–600	[2, 43]
μ	Rate at which children leave age group below 12 or die	$2\times {10}^{-4}-5\times {10}^{-4}$	[2, 23, 43]
β	Usual transmission coefficient	0.1–1	[21, 33]
${\alpha }_1$	Relative contact rate of symptomatic cases in the presence of constant behaviour change	0.09–0.5	[34, 43]
${\alpha }_2$	Relative contact rate of quarantined individuals in the presence of constant behaviour change	0.01–0.4	Assumed
${\sigma }^{-1}$	incubation period	2–5	[38]
κ	Rate at which symptomatic cases are quarantined	0.1–1	[27]
ρ	Proportion of infections being symptomatic	0.4–0.6	[44]
${\gamma }_1$	Recovery rate of symptomatic infections	0.08–0.2	[21, 33, 38]
${\gamma }_2$	Recovery rate of asymptomatic infections	0.08–0.2	[21, 43]
${\gamma }_3$	Recovery rate of quarantined individuals	0.08–0.3	[21, 27]

Proof.

The vector field generated by the right side of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) is continuously differentiable in Ω, so it is Lipschitz continuous. Thus, given any nonnegative initial condition, there exists a unique solution for all $t\ge 0$. If S = 0, then $S^{\prime }=\mathrm{\Lambda }>0$. It follows from Proposition B.7 in Smith and Waltman [36] that $S(t)\ge 0$. The nonnegativity of other state variables can be similarly proved. If $S+E+I_s+I_a+Q=\mathrm{\Lambda }/\mu$, then $S^{\prime }+E^{\prime }+I_s^{\prime }+I_a^{\prime }+Q^{\prime }=-({\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q)\le 0$. This means that Ω is positively invariant with respect to system (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ).

3. Mathematical analysis

In this section, we use the next generation matrix method [11] to define the basic reproduction number of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) and then establish the threshold dynamic result for the model.

3.1. Basic reproduction number

Setting the right-hand side of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) to zero gives a unique disease-free equilibrium ${\mathit{E}}_0=(\mathrm{\Lambda }/\mu ,0,0,0,0).$Using the recipe of van den Driessche and Watmough [39], we get $\mathcal{F}=\left(\begin{array}{c}\lambda (t)S\\ 0\\ 0\\ 0\end{array}\right)\text{and}\mathcal{V}=\left(\begin{array}{c}(\sigma +\mu )E\\ -\rho \sigma E+(\kappa +{\gamma }_1+\mu )I_s\\ -(1-\rho )\sigma E+({\gamma }_2+\mu )I_a\\ -\kappa I_s+({\gamma }_3+\mu )Q\end{array}\right).$Linearizing $\mathcal{F}$ and $\mathcal{V}$ at ${\mathit{E}}_0$ gives the incidence and transition matrices $\begin{array}{rl}F& =\left(\begin{array}{cccc}0& {\alpha }_1\beta & \beta & {\alpha }_2\beta \\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right)\text{and}\\ V& =\left(\begin{array}{cccc}\sigma +\mu & 0& 0& 0\\ -\rho \sigma & \kappa +{\gamma }_1+\mu & 0& 0\\ -(1-\rho )\sigma & 0& {\gamma }_2+\mu & 0\\ 0& -\kappa & 0& {\gamma }_3+\mu \end{array}\right),\end{array}$respectively. The basic reproduction number of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) is defined as the spectral radius of the next generation matrix $F{V}^{-1}$, i.e. ${\mathcal{R}}_0=\rho (F{V}^{-1})={\mathcal{R}}_1+{\mathcal{R}}_2+{\mathcal{R}}_3,$where $\begin{array}{rl}{\mathcal{R}}_1& =\frac{\sigma }{\sigma +\mu }\cdot \rho \cdot \beta {\alpha }_1\cdot \frac{1}{\kappa +{\gamma }_1+\mu },\\ {\mathcal{R}}_2& =\frac{\sigma }{\sigma +\mu }\cdot (1-\rho )\cdot \beta \cdot \frac{1}{{\gamma }_2+\mu },\\ {\mathcal{R}}_3& =\frac{\sigma }{\sigma +\mu }\cdot \rho \cdot \frac{\kappa }{\kappa +{\gamma }_1+\mu }\cdot \beta {\alpha }_2\cdot \frac{1}{{\gamma }_3+\mu }\end{array}$represent the average number of secondary cases produced by one infected individual in the person's symptomatic, asymptomatic and quarantined stages, respectively. In ${\mathcal{R}}_1$, the ratio $\sigma /(\sigma +\mu )$ is the probability that an exposed individual will survive the incubation period and become infectious, ρ is the fraction of exposed individuals who will develop symptoms, $\beta {\alpha }_1$ is the transmission coefficient for symptomatically infected individuals, and $1/(\kappa +{\gamma }_1+\mu )$ is the average infectious period of a symptomatic individual. In ${\mathcal{R}}_2$, the term $1-\rho$ is the fraction of exposed individuals who will be symptomless, β is the transmission coefficient for the asymptomatically infected individuals, and $1/({\gamma }_2+\mu )$ is the average infectious period of an asymptomatic individual. In ${\mathcal{R}}_3$, the ratio $\kappa /(\kappa +{\gamma }_1+\mu )$ is the probability that a symptomatic individual will be quarantined, $\beta {\alpha }_2$ is the transmission coefficient for quarantined individuals, and $1/({\gamma }_3+\mu )$ is the average infectious period of a quarantined individual. Clearly, only behaviour responses from symptomatically infected and quarantined people affect the basic reproduction number.

3.2. Threshold dynamics

In what follows, we study the disease dynamics of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) by Lyapunov method.

Theorem 3.1

For model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ), if ${\mathcal{R}}_0\le 1$, then the disease-free equilibrium ${\mathit{E}}_0$ is globally asymptotically stable; if ${\mathcal{R}}_0>1$, then the disease-free equilibrium ${\mathit{E}}_0$ is unstable.

Proof.

Following Theorem 2 in van den Driessche and Watmough [39], the local asymptotic stability of the disease-free equilibrium ${\mathit{E}}_0$ can be immediately obtained. Thus, it suffices to show the global attractivity of ${\mathit{E}}_0$ in Ω as ${\mathcal{R}}_0\le 1$. We construct a Lyapunov function $L=E+h_1{I}_s+h_2{I}_a+h_{3}Q,$where constants $h_1,h_2$ and $h_3$ are yet to be determined. Differentiating L along (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) gives $\begin{array}{rl}{L}^{\prime }& =E^{\prime }+h_1{I}_s^{\prime }+h_2{I}_a^{\prime }+h_3{Q}^{\prime }\\ & =\beta {\displaystyle \frac{{\alpha }_1{I}_s+f(I_s,Q)I_a+{\alpha }_{2}Q}{N^{\ast }}}S-(\sigma +\mu )E+h_1{I}_s^{\prime }+h_2{I}_a^{\prime }+h_3{Q}^{\prime }.\end{array}$The facts $S\le N^{\ast }$ and $f(I_s,Q)\le 1$ indicate that $\begin{array}{rl}{L}^{\prime }& \le \beta ({\alpha }_1{I}_s+I_a+{\alpha }_{2}Q)-(\sigma +\mu )E+h_1(\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s)\\ & +h_2((1-\rho )\sigma E-({\gamma }_2+\mu )I_a)+h_3(\kappa I_s-({\gamma }_3+\mu )Q).\end{array}$Combining like terms yields (3) $\begin{array}{rl}{L}^{\prime }& \le (h_1\rho \sigma +h_2(1-\rho )\sigma -(\sigma +\mu ))E+(\beta {\alpha }_1-h_1(\kappa +{\gamma }_1+\mu )+h_3\kappa )I_s\\ & +(\beta -h_2({\gamma }_2+\mu ))I_a+(\beta {\alpha }_2-h_3({\gamma }_3+\mu ))Q.\end{array}$(3) To find appropriate constants $h_1,h_2$ and $h_3$ so that the last three terms on the right-hand side of (3(3) $\begin{array}{rl}{L}^{\prime }& \le (h_1\rho \sigma +h_2(1-\rho )\sigma -(\sigma +\mu ))E+(\beta {\alpha }_1-h_1(\kappa +{\gamma }_1+\mu )+h_3\kappa )I_s\\ & +(\beta -h_2({\gamma }_2+\mu ))I_a+(\beta {\alpha }_2-h_3({\gamma }_3+\mu ))Q.\end{array}$(3) ) equal zero, we set $\begin{array}{rl}& \beta {\alpha }_1-h_1(\kappa +{\gamma }_1+\mu )+h_3\kappa =0,\\ & \beta -h_2({\gamma }_2+\mu )=0,\\ & \beta {\alpha }_2-h_3({\gamma }_3+\mu )=0.\end{array}$Solving the above gives $(h_1,h_2,h_3)=({\displaystyle \frac{\beta {\alpha }_1}{\kappa +{\gamma }_1+\mu }}+{\displaystyle \frac{\kappa \beta {\alpha }_2}{(\kappa +{\gamma }_1+\mu )({\gamma }_3+\mu )}},{\displaystyle \frac{\beta }{{\gamma }_2+\mu }},{\displaystyle \frac{\beta {\alpha }_2}{{\gamma }_3+\mu }}).$Substituting the obtained $h_1,h_2$ and $h_3$ into (3(3) $\begin{array}{rl}{L}^{\prime }& \le (h_1\rho \sigma +h_2(1-\rho )\sigma -(\sigma +\mu ))E+(\beta {\alpha }_1-h_1(\kappa +{\gamma }_1+\mu )+h_3\kappa )I_s\\ & +(\beta -h_2({\gamma }_2+\mu ))I_a+(\beta {\alpha }_2-h_3({\gamma }_3+\mu ))Q.\end{array}$(3) ), we get $\begin{array}{rl}{L}^{\prime }& \le ({\displaystyle \frac{\rho \sigma \beta {\alpha }_1}{\kappa +{\gamma }_1+\mu }}+{\displaystyle \frac{(1-\rho )\sigma \beta }{{\gamma }_2+\mu }}+\frac{\rho \sigma \kappa \beta {\alpha }_2}{(\kappa +{\gamma }_1+\mu )({\gamma }_3+\mu )}-(\sigma +\mu ))E\\ & =({\mathcal{R}}_0-1)(\sigma +\mu )E.\end{array}$If ${\mathcal{R}}_0<1$, then $L^{\prime }\le 0$ and $L^{\prime }=0$ if and only if E = 0. The second and third equations of (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) give $I_s^{\prime }=-(\kappa +{\gamma }_1+\mu )I_s\text{and}{I}_a^{\prime }=-({\gamma }_2+\mu )I_a,$which imply $I_s=I_a=0$. Thus, the last equation of (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) gives $Q^{\prime }=-({\gamma }_3+\mu )Q$ and hence Q = 0. The first equation of (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) becomes $S^{\prime }=\mathrm{\Lambda }-\mu S$ and thus $S=\mathrm{\Lambda }/\mu$. If ${\mathcal{R}}_0=1$, then the equality $L^{\prime }=0$ implies that $\beta {\displaystyle \frac{{\alpha }_1{I}_s+f(I_s,Q)I_a+{\alpha }_{2}Q}{N^{\ast }}}S=\beta ({\alpha }_1{I}_s+I_a+{\alpha }_{2}Q).$Then $S=N^{\ast }$ and $I_a=0$ or $I_s=Q=0$, or $I_s=I_a=Q=0$ holds and we can proceed as before. In summary, if ${\mathcal{R}}_0\le 1$, then the largest compact invariant subset of $\{(S,E,I_s,I_a,Q)\in \mathrm{\Omega }|L^{\prime }=0\}$ is the singleton $\{{\mathit{E}}_0\}$. By the LaSalle's invariance principle, the disease-free equilibrium ${\mathit{E}}_0$ is globally asymptotically stable as ${\mathcal{R}}_0\le 1$.

Since the matrix $V^{-1}F$ associated to model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) is reducible, Theorem 2.2 in Shuai and van den Driessche [35] does not work directly. However, by a suitable split of the matrix F−V (see Subsection 3.3 in Gao and Cao [14]), we can still use a similar approach to construct an implicit Lyapunov function. Next we consider the disease transmission under ${\mathcal{R}}_0>1$.

Theorem 3.2

Assume that ${\mathcal{R}}_0>1$, then system (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) is uniformly persistent, i.e. there is a constant $\epsilon >0$ such that each solution ${\phi }_t({\mathit{x}}_0)=(S(t),E(t),I_s(t),I_a(t),Q(t))$ with the initial value ${\mathit{x}}_0=(S(0),E(0),I_s(0),I_a(0),Q(0))\in \mathrm{\Omega }$ satisfies $\underset{t\to \mathrm{\infty }}{lim inf}(E(t),I_s(t),I_a(t),Q(t))\gg (\epsilon ,\epsilon ,\epsilon ,\epsilon ),$where $E(0)+I_s(0)+I_a(0)+Q(0)>0$.

Proof.

We use Theorem 4.6 in Thieme [37] to prove the persistence. Denote $\begin{array}{rl}{\mathrm{\Omega }}_0& =\{(S,E,I_s,I_a,Q)\in \mathrm{\Omega }|E>0,I_s>0,I_a>0,Q>0\},\\ \mathrm{\partial }{\mathrm{\Omega }}_0& =\mathrm{\Omega }\mathrm{\setminus }{\mathrm{\Omega }}_0=\{(S,E,I_s,I_a,Q)\in \mathrm{\Omega }|E=0\text{or}{I}_s=0\text{or}{I}_a=0\text{or}Q=0\}.\end{array}$Obviously, Ω and ${\mathrm{\Omega }}_0$ are positively invariant and $\mathrm{\partial }{\mathrm{\Omega }}_0$ is a relatively closed subset of Ω. The point dissipativity of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) can be seen from Theorem 2.1. Set $\begin{array}{rl}{M}_{\mathrm{\partial }}& =\{{\mathit{x}}_0\in \mathrm{\partial }{\mathrm{\Omega }}_0|{\phi }_t({\mathit{x}}_0)\in \mathrm{\partial }{\mathrm{\Omega }}_0,\mathrm{\forall }t\ge 0\},\\ D& =\{(S,0,0,0,0)\in \mathrm{\Omega }|S\ge 0\}.\end{array}$We claim that $M_{\mathrm{\partial }}=D$. Note that $D\subseteq M_{\mathrm{\partial }}$, so we only need to prove $M_{\mathrm{\partial }}\subseteq D$. If ${\mathit{x}}_0\in \mathrm{\partial }{\mathrm{\Omega }}_0\mathrm{\setminus }D$, then $E(0)+I_s(0)+I_a(0)+Q(0)>0$. It follows from the irreducibility of F−V that ${\phi }_t({\mathit{x}}_0)\in {\mathrm{\Omega }}_0$ for $\mathrm{\forall }t>0$. Therefore, we have ${\phi }_t({\mathit{x}}_0)\notin \mathrm{\partial }{\mathrm{\Omega }}_0$ and ${\mathit{x}}_0\notin M_{\mathrm{\partial }}$ which implies $M_{\mathrm{\partial }}=D$.

The only equilibrium contained in $M_{\mathrm{\partial }}=D$ is ${\mathit{E}}_0$. So $\bigcup _{{\mathit{x}}_0\in M_{\mathrm{\partial }}}\omega ({\mathit{x}}_0)=\{{\mathit{E}}_0\}$. Thus, $\{{\mathit{E}}_0\}$ is a compact and isolated invariant set in $M_{\mathrm{\partial }}$. Let $W^s({\mathit{E}}_0)$ be the stable manifold of ${\mathit{E}}_0$. It remains to prove that $W^s({\mathit{E}}_0)\bigcap {\mathrm{\Omega }}_0=\mathrm{\varnothing }$ when ${\mathcal{R}}_0>1$. Suppose not, then there exists ${\mathit{x}}_0\in {\mathrm{\Omega }}_0$ such that ${\phi }_t({\mathit{x}}_0)\to {\mathit{E}}_0$ as $t\to +\mathrm{\infty }$. Define $M(\xi )=F-V-\xi K=\left(\begin{array}{cccc}-(\sigma +\mu )& (1-2\xi ){\alpha }_1\beta & (1-2\xi )\beta & (1-2\xi ){\alpha }_2\beta \\ \rho \sigma & -(\kappa +{\gamma }_1+\mu )& 0& 0\\ (1-\rho )\sigma & 0& -({\gamma }_2+\mu )& 0\\ 0& \kappa & 0& -({\gamma }_3+\mu )\end{array}\right),$where $K=\left(\begin{array}{cccc}0& 2{\alpha }_1\beta & 2\beta & 2{\alpha }_2\beta \\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right).$Using Theorem 2 in van den Driessche and Watmough [39], we know the spectral bound $s(M(0))=s(F-V)>0$ if and only if ${\mathcal{R}}_0=\rho (F{V}^{-1})>1$. Since $s(M(\xi ))$ is continuous in ξ, there exists a small ${\xi }_1$ such that $s(M(\xi ))>0$ for $\xi \in [0,{\xi }_1]$.

It follows from $\underset{t\to \mathrm{\infty }}{lim}{\phi }_t({\mathit{x}}_0)={\mathit{E}}_0$ that, for some T>0, we have ${\Vert {\phi }_t({\mathit{x}}_0)-{\mathit{E}}_0\Vert }_2\le \tau ,\mathrm{\forall }t>T,$where $\Vert \cdot {\Vert }_2$ is the Euclidean norm and τ is small enough such that $\frac{S(t)}{N^{\ast }}\ge 1-{\xi }_1\text{and}f(I_s(t),Q(t))\ge 1-{\xi }_1.$Therefore, when t>T, we have $\begin{array}{rl}\frac{\mathrm{d}E}{\mathrm{d}t}& \ge ({\alpha }_1\beta I_s+(1-{\xi }_1)\beta I_a+{\alpha }_2\beta Q)(1-{\xi }_1)-(\sigma +\mu )E\\ & \ge (1-2{\xi }_1)({\alpha }_1\beta I_s+\beta I_a+{\alpha }_2\beta Q)-(\sigma +\mu )E,\end{array}$which implies $\begin{array}{rl}\frac{\mathrm{d}E}{\mathrm{d}t}& \ge (1-2{\xi }_1)({\alpha }_1\beta I_s+\beta I_a+{\alpha }_2\beta Q)-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q.\end{array}$Consider an auxiliary system (4) ${\mathit{u}}^{\prime }=M({\xi }_1)\mathit{u},$(4) where $\mathit{u}=(u_1,u_2,u_3,u_4{)}^{\mathrm{T}}=(E,I_s,I_a,Q{)}^{\mathrm{T}}$. Since $M({\xi }_1)$ is essentially nonnegative and irreducible, by the Perron–Frobenius theorem, it has a positive eigenvector associated to $s(M({\xi }_1))>0$. Thus, any positive solution of system (4(4) ${\mathit{u}}^{\prime }=M({\xi }_1)\mathit{u},$(4) ) satisfies $u_i(t)\to +\mathrm{\infty }$ as $t\to +\mathrm{\infty }$, i = 1, 2, 3, 4. Applying the comparison principle [36], we get $\underset{t\to +\mathrm{\infty }}{lim}(E(t),I_s(t),I_a(t),Q(t))=(+\mathrm{\infty },+\mathrm{\infty },+\mathrm{\infty },+\mathrm{\infty }).$This gives a contradiction. So $W^s({\mathit{E}}_0)\bigcap {\mathrm{\Omega }}_0=\mathrm{\varnothing }$. By Theorem 4.6 in Thieme [37], the system is uniformly persistent if ${\mathcal{R}}_0>1$.

Thus, ${\mathcal{R}}_0$ is a threshold quantity that determines disease extinction and persistence. By Theorem 2.4 in Zhao [47], the uniform persistence of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) implies the existence of at least one endemic equilibrium. Furthermore, we can prove the uniqueness of endemic equilibrium in a direct way.

Theorem 3.3

For model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ), there exists exactly one endemic equilibrium if and only if ${\mathcal{R}}_0>1$.

Proof.

By Theorem 3.1, it suffices to show the uniqueness of endemic equilibrium. Denote the endemic equilibrium of (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) by ${\mathit{E}}^{\ast }=(S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })$ which satisfies (5a) $\begin{array}{rl}& \mathrm{\Lambda }-\lambda (S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })S^{\ast }-\mu S^{\ast }=0,\end{array}$(5a) (5b) $\begin{array}{rl}& \lambda (S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })S^{\ast }-(\sigma +\mu )E^{\ast }=0,\end{array}$(5b) (5c) $\begin{array}{rl}& \rho \sigma E^{\ast }-(\kappa +{\gamma }_1+\mu )I_s^{\ast }=0,\end{array}$(5c) (5d) $\begin{array}{rl}& (1-\rho )\sigma E^{\ast }-({\gamma }_2+\mu )I_a^{\ast }=0,\end{array}$(5d) (5e) $\begin{array}{rl}& \kappa I_s^{\ast }-({\gamma }_3+\mu )Q^{\ast }=0.\end{array}$(5e) Solving $E^{\ast },I_a^{\ast }$ and $Q^{\ast }$ from (5c(5c) $\begin{array}{rl}& \rho \sigma E^{\ast }-(\kappa +{\gamma }_1+\mu )I_s^{\ast }=0,\end{array}$(5c) )–(5e(5e) $\begin{array}{rl}& \kappa I_s^{\ast }-({\gamma }_3+\mu )Q^{\ast }=0.\end{array}$(5e) ) gives (6) $\begin{array}{rl}{E}^{\ast }& =\frac{\kappa +{\gamma }_1+\mu }{\rho \sigma }{I}_s^{\ast },\\ I_a^{\ast }& =\frac{(1-\rho )(\kappa +{\gamma }_1+\mu )}{\rho ({\gamma }_2+\mu )}{I}_s^{\ast },\\ Q^{\ast }& =\frac{\kappa }{{\gamma }_3+\mu }{I}_s^{\ast }.\end{array}$(6) Substituting (6(6) $\begin{array}{rl}{E}^{\ast }& =\frac{\kappa +{\gamma }_1+\mu }{\rho \sigma }{I}_s^{\ast },\\ I_a^{\ast }& =\frac{(1-\rho )(\kappa +{\gamma }_1+\mu )}{\rho ({\gamma }_2+\mu )}{I}_s^{\ast },\\ Q^{\ast }& =\frac{\kappa }{{\gamma }_3+\mu }{I}_s^{\ast }.\end{array}$(6) ) into (5b(5b) $\begin{array}{rl}& \lambda (S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })S^{\ast }-(\sigma +\mu )E^{\ast }=0,\end{array}$(5b) ) yields ${\displaystyle \frac{\beta S^{\ast }{I}_s^{\ast }}{N^{\ast }}}({\alpha }_1+f(I_s^{\ast },Q^{\ast }){\displaystyle \frac{(1-\rho )(\kappa +{\gamma }_1+\mu )}{\rho ({\gamma }_2+\mu )}}+{\displaystyle \frac{{\alpha }_2\kappa }{{\gamma }_3+\mu }})=(\sigma +\mu ){\displaystyle \frac{\kappa +{\gamma }_1+\mu }{\rho \sigma }}{I}_s^{\ast },$which can be rewritten as ${\displaystyle \frac{S^{\ast }}{N^{\ast }}}\cdot {\displaystyle \frac{(\sigma +\mu )(\kappa +{\gamma }_1+\mu )}{\rho \sigma }}\cdot ({\mathcal{R}}_1+f(I_s^{\ast },Q^{\ast }){\mathcal{R}}_2+{\mathcal{R}}_3)={\displaystyle \frac{(\sigma +\mu )(\kappa +{\gamma }_1+\mu )}{\rho \sigma }}.$Thus, we have (7) $S^{\ast }={\displaystyle \frac{N^{\ast }}{{\mathcal{R}}_1+f(I_s^{\ast },Q^{\ast }){\mathcal{R}}_2+{\mathcal{R}}_3}}.$(7) The sum of (5a(5a) $\begin{array}{rl}& \mathrm{\Lambda }-\lambda (S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })S^{\ast }-\mu S^{\ast }=0,\end{array}$(5a) ) and (5b(5b) $\begin{array}{rl}& \lambda (S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })S^{\ast }-(\sigma +\mu )E^{\ast }=0,\end{array}$(5b) ) is (8) $\mathrm{\Lambda }-\mu S^{\ast }=(\sigma +\mu )E^{\ast }.$(8) Substituting (6(6) $\begin{array}{rl}{E}^{\ast }& =\frac{\kappa +{\gamma }_1+\mu }{\rho \sigma }{I}_s^{\ast },\\ I_a^{\ast }& =\frac{(1-\rho )(\kappa +{\gamma }_1+\mu )}{\rho ({\gamma }_2+\mu )}{I}_s^{\ast },\\ Q^{\ast }& =\frac{\kappa }{{\gamma }_3+\mu }{I}_s^{\ast }.\end{array}$(6) ) and (7(7) $S^{\ast }={\displaystyle \frac{N^{\ast }}{{\mathcal{R}}_1+f(I_s^{\ast },Q^{\ast }){\mathcal{R}}_2+{\mathcal{R}}_3}}.$(7) ) into (8(8) $\mathrm{\Lambda }-\mu S^{\ast }=(\sigma +\mu )E^{\ast }.$(8) ) yields $\begin{array}{rl}{I}_s^{\ast }& =\frac{\mathrm{\Lambda }-\mu S^{\ast }}{\sigma +\mu }\cdot \frac{\rho \sigma }{\kappa +{\gamma }_1+\mu }\\ & =(\mathrm{\Lambda }-\frac{\mu N^{\ast }}{{\mathcal{R}}_1+f(I_s^{\ast },Q^{\ast }){\mathcal{R}}_2+{\mathcal{R}}_3})\frac{\rho \sigma }{(\sigma +\mu )(\kappa +{\gamma }_1+\mu )}.\end{array}$It follows that $f(I_s^{\ast },Q^{\ast })=(\frac{a\mathrm{\Lambda }}{a\mathrm{\Lambda }-I_s^{\ast }}-{\mathcal{R}}_1-{\mathcal{R}}_3)\frac{1}{{\mathcal{R}}_2},$where $a={\displaystyle \frac{\rho \sigma }{(\sigma +\mu )(\kappa +{\gamma }_1+\mu )}}.$Denote $G(I_s^{\ast })=f(I_s^{\ast },{\displaystyle \frac{\kappa }{{\gamma }_3+\mu }}{I}_s^{\ast })-({\displaystyle \frac{a\mathrm{\Lambda }}{a\mathrm{\Lambda }-I_s^{\ast }}}-{\mathcal{R}}_1-{\mathcal{R}}_3)\frac{1}{{\mathcal{R}}_2},I_s^{\ast }\in [0,a\mathrm{\Lambda }).$Differentiating $G(I_s^{\ast })$ with respect to $I_s^{\ast }$, we find that $G(I_s^{\ast })$ is strictly decreasing in $I_s^{\ast }\in [0,a\mathrm{\Lambda })$. Note that $G(0)=1-(1-{\mathcal{R}}_1-{\mathcal{R}}_3)\frac{1}{{\mathcal{R}}_2}={\displaystyle \frac{{\mathcal{R}}_1+{\mathcal{R}}_2+{\mathcal{R}}_3-1}{{\mathcal{R}}_2}}={\displaystyle \frac{{\mathcal{R}}_0-1}{{\mathcal{R}}_2}}>0\text{ as }{\mathcal{R}}_0>1$and $\underset{I_s^{\ast }\to a\mathrm{\Lambda }-}{lim}G(I_s^{\ast })=-\mathrm{\infty }\text{ and }a\mathrm{\Lambda }<N^{\ast }=\frac{\mathrm{\Lambda }}{\mu }.$Thus, the equation $G(I_s^{\ast })=0$ has a unique positive root when ${\mathcal{R}}_0>1$.

Remark 3.4

For simple endemic models like SIS, SIAR and SEIRS, the fraction of people being susceptible at the endemic equilibrium (when it exists) usually equals the reciprocal of the basic reproduction number, i.e. $S^{\ast }/N^{\ast }=1/{\mathcal{R}}_0$. However, for model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ), it follows from (7(7) $S^{\ast }={\displaystyle \frac{N^{\ast }}{{\mathcal{R}}_1+f(I_s^{\ast },Q^{\ast }){\mathcal{R}}_2+{\mathcal{R}}_3}}.$(7) ) that $\frac{S^{\ast }}{N^{\ast }}={\displaystyle \frac{1}{{\mathcal{R}}_1+f(I_s^{\ast },Q^{\ast }){\mathcal{R}}_2+{\mathcal{R}}_3}}\ge {\displaystyle \frac{1}{{\mathcal{R}}_1+{\mathcal{R}}_2+{\mathcal{R}}_3}}={\displaystyle \frac{1}{{\mathcal{R}}_0}}.$Biologically speaking, behaviour change reduces the number of infections. Indeed, (6(6) $\begin{array}{rl}{E}^{\ast }& =\frac{\kappa +{\gamma }_1+\mu }{\rho \sigma }{I}_s^{\ast },\\ I_a^{\ast }& =\frac{(1-\rho )(\kappa +{\gamma }_1+\mu )}{\rho ({\gamma }_2+\mu )}{I}_s^{\ast },\\ Q^{\ast }& =\frac{\kappa }{{\gamma }_3+\mu }{I}_s^{\ast }.\end{array}$(6) ) indicates that the population size of each nonsusceptible compartment at the endemic equilibrium ${\mathit{E}}^{\ast }$ becomes smaller when behaviour change occurs. Moreover, it weakens the relative contribution of asymptomatic class in disease transmission. In fact, the ratio of the forces of infections attributed to $I_s,I_a$ and Q at ${\mathit{E}}^{\ast }$ is $\begin{array}{rl}& \frac{{\alpha }_1\beta I_s^{\ast }}{N^{\ast }}:\frac{f(I_s^{\ast },Q^{\ast })\beta I_a^{\ast }}{N^{\ast }}:\frac{{\alpha }_2\beta Q^{\ast }}{N^{\ast }}={\alpha }_1{I}_s^{\ast }:f(I_s^{\ast },Q^{\ast })I_a^{\ast }:{\alpha }_2{Q}^{\ast }\\ & ={\alpha }_1\frac{\rho \sigma }{\kappa +{\gamma }_1+\mu }{E}^{\ast }:f(I_s^{\ast },Q^{\ast })\frac{(1-\rho )\sigma }{{\gamma }_2+\mu }{E}^{\ast }:{\alpha }_2\frac{\kappa }{{\gamma }_3+\mu }\cdot \frac{\rho \sigma }{\kappa +{\gamma }_1+\mu }{E}^{\ast }\\ & \le \frac{{\alpha }_1\rho }{\kappa +{\gamma }_1+\mu }:\frac{1-\rho }{{\gamma }_2+\mu }:\frac{\kappa }{\kappa +{\gamma }_1+\mu }\cdot \frac{{\alpha }_2\rho }{{\gamma }_3+\mu }={\mathcal{R}}_1:{\mathcal{R}}_2:{\mathcal{R}}_3,\end{array}$where the second equality is due to the equilibrium Equation (5). Interestingly, in the absence of behaviour change, the relative contributions of infected states $I_s,I_a$ and Q to new infections at the endemic equilibrium are the same as their contributions to the basic reproduction number.

Before ending this section, we give a sufficient condition under which the endemic equilibrium is globally attractive. In particular, it is satisfied when $f(I_s,Q)$ is constant. The proof is postponed to Appendix A.

Theorem 3.5

Suppose that $(f(I_s,Q)-f(I_s^{\ast },Q^{\ast }))(f(I_s,Q)I_a-f(I_s^{\ast },Q^{\ast })I_a^{\ast })\le 0$holds for $I_s,I_a,Q>0$ and $I_s+I_a+Q\le N^{\ast }=\mathrm{\Lambda }/\mu$. If ${\mathcal{R}}_0>1$, then the unique endemic equilibrium ${\mathit{E}}^{\ast }$ of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) is globally asymptotically stable in Ω minus the disease-free space.

4. Numerical analysis

In this part, we first fit the proposed epidemiological model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q,R)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q,R)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& ={\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q-\mu R,\end{array}$(1) ) to weekly HFMD reported case data in Shanghai, China. Then, we carry out some sensitivity analysis and numerical simulations to compare the effectiveness of different intervention strategies and explore the role of behaviour change on infection control.

4.1. Date fitting

The first HFMD case in mainland China was emerged in Shanghai in 1981. It has been mandatory to report cases of HFMD to Shanghai Municipal Center for Disease Control and Prevention since 2005, three years earlier than the country. Shanghai has established a city-wide communicable disease surveillance system that consists of medical institutions, district and municipal centres for disease control and prevention. Every year tens of thousands of cases are reported in Shanghai and almost all of them are admitted to three specialized hospitals. The usual peak season in Shanghai is from April to July and a smaller peak may occur from September to November. Although Shanghai is one of the richest cities in China, its annual HFMD incidence rate is significantly higher than the national average partially due to meteorological factors and migrant population. For example, the incidence rates of Shanghai and China in 2016 are $237/100,000$ and $178/100,000$, respectively. Schools and kindergartens implement emergency response plan for childhood infectious diseases and class or school closure is common when a cluster of cases are identified in a short time period.

By using the least-squares method, we fit model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q,R)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q,R)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& ={\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q-\mu R,\end{array}$(1) ) to the weekly HFMD case data of Shanghai in 2019 (Zhao et al. [49]). According to the 2019 Shanghai Statistical Yearbook [2], we set $\mathrm{\Lambda }=463$, $\mu =2.3\times {10}^{-4}$ and $N(0)=2.01\times {10}^6$. To reduce uncertainty, we choose some reasonable parameter values from literature as follows ${\alpha }_1=0.5,{\alpha }_2=0.1,\rho =0.4,\sigma =0.25,\kappa =0.5,\text{and}{\gamma }_3=0.25,$where the time unit is one day. We assume that ${\gamma }_1={\gamma }_2=\gamma$, indicating that individuals either symptomatic or asymptomatic will take ${\gamma }^{-1}$ days to recover. Since it takes time to put people into quarantine, it is reasonably believed that ${\gamma }^{-1}>{\gamma }_3^{-1}$, i.e. $\gamma <{\gamma }_3$. The behaviour change function is assumed to take the form (9) $f(I_s,Q)={\displaystyle \frac{1}{1+m_1{I}_s+m_{2}Q}},$(9) where $m_1>0$ and $m_2>0$ measure the influence of newly reported cases and recently cumulative cases on the contact rate of ‘healthy people’, respectively. The initial conditions of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q,R)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q,R)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& ={\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q-\mu R,\end{array}$(1) ) are set as $(S(0),E(0),I_s(0),I_a(0),Q(0),R(0))=(1.01\times {10}^6,635,91,682,350,{10}^6).$Denote the given data set by $\{(t_1,Y_1),\dots ,(t_n,Y_n)\}$ where $Y_i$ is the newly reported number of cases in week i. The simulated weekly cases are obtained as $Z_i={\int }_{t_{i-1}}^{t_i}\kappa I_s(t)\mathrm{d}t,$and the sum-of-squares error is defined as: $\text{SSE}=\sum _{i=1}^n(Y_i-Z_i{)}^2.$We use the open-source R programming language to calibrate the model and estimate the disease transmission coefficient, β, behaviour change parameters, $m_1$ and $m_2$, and recovery rate of infection, γ. Specifically, based on the Levenberg–Marquardt algorithm, we use the nls.lm function from the R package minpack.lm to perform the least-squares fitting. Note that in our model calibration, the transmission coefficient, β, is considered as a time-dependent cubic B-spline function instead of a constant to better capture the seasonality in reported data. The fitted parameter values are $m_1=3.8\times {10}^{-4}$, $m_2=7.2\times {10}^{-4}$, and $\gamma =0.14$. Figure 3 illustrates the simulated average weekly cases versus time, and the grey shaded area gives the 95% confidence interval for the number of simulated cases per week, which matches the changing pattern of the reported weekly cases. Hence our model provides a good fit for the time series of the weekly reported HFMD cases. Moreover, the time evolution of the basic reproduction number is estimated and shown in Figure 4, which indicates that the asymptotically infected individuals contributes the most to the transmission. This is mainly attributed to the rapid and strict quarantine of symptomatic cases. The average basic reproduction number and transmission coefficient are ${\overline{\mathcal{R}}}_0=2.41$ and $\overline{\beta }=0.52$, respectively, which indicate that HFMD cannot be eradicated under the current control strategy.

Figure 3. Comparison of reported and simulated weekly HFMD cases of Shanghai in 2019. The red triangles are reported weekly cases, the blue dots are simulated weekly cases, and the grey shaded area gives the 95% confidence interval for the number of simulated cases per week.

Figure 4. Time evolution of ${\mathcal{R}}_0={\mathcal{R}}_1+{\mathcal{R}}_2+{\mathcal{R}}_3$, where ${\mathcal{R}}_1,{\mathcal{R}}_2$ and ${\mathcal{R}}_3$ denote the secondary infections produced by one infected individual during symptomatic, asymptomatic and quarantined stages, respectively.

4.2. Sensitivity analysis

Since the threshold dynamics of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) are wholly governed by the basic reproduction number, reducing ${\mathcal{R}}_0$ to be less than one is desirable for disease eradication. Thus, it is necessary to examine how ${\mathcal{R}}_0$ varies with model parameters. Direct calculations find that ${\mathcal{R}}_0$ is monotone increasing with respect to $\beta ,{\alpha }_1,{\alpha }_2$ and σ, decreasing in ${\gamma }_1,{\gamma }_2,{\gamma }_3$ and μ, and independent of Λ and f. Mathematically, the monotonicity of ${\mathcal{R}}_0$ on ρ and κ may change with parameter setting by noting that $\mathrm{sgn}\left(\frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }\rho }\right)=\mathrm{sgn}(\frac{{\alpha }_1}{\kappa +{\gamma }_1+\mu }-\frac{1}{{\gamma }_2+\mu }+\frac{\kappa }{\kappa +{\gamma }_1+\mu }\cdot \frac{{\alpha }_2}{{\gamma }_3+\mu })$and $\mathrm{sgn}\left(\frac{\mathrm{\partial }{\mathcal{R}}_0}{\mathrm{\partial }\kappa }\right)=\mathrm{sgn}(\frac{{\alpha }_2}{{\gamma }_3+\mu }-\frac{{\alpha }_1}{{\gamma }_1+\mu }).$In the real world, the highly possible relations ${\alpha }_2<{\alpha }_1<1$, ${\gamma }_1<{\gamma }_3$ and ${\gamma }_2<\kappa +{\gamma }_1$ imply that ${\mathcal{R}}_0$ is probably decreasing in ρ and κ.

To see which way is the best in lowering disease persistence, we have to compare the sensitivity indices of ${\mathcal{R}}_0$ to parameter variation. The normalized forward sensitivity index is widely used and it is defined as the ratio of the relative change of the output to the relative change of the parameter [3]. Specifically, let u be a variable that depends differentially on a parameter s. Then the sensitivity index (SI) of u in terms of s is expressed as ${\mathrm{{\rm Y}}}_s^u=\underset{\mathrm{\Delta }s\to 0}{lim}{\displaystyle \frac{\mathrm{\Delta }u}{u}}/{\displaystyle \frac{\mathrm{\Delta }s}{s}}=\underset{\mathrm{\Delta }s\to 0}{lim}{\displaystyle \frac{\mathrm{\Delta }u}{\mathrm{\Delta }s}}\cdot {\displaystyle \frac{s}{u}}=\frac{\mathrm{\partial }u}{\mathrm{\partial }s}\cdot \frac{s}{u}.$The complex structure of ${\mathcal{R}}_0$ for model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q,R)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q,R)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& ={\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q-\mu R,\end{array}$(1) ) makes it difficult to analytically compare the sensitivity indices of different parameters. Therefore, we select a typical set of parameter values based on the above fitting result: (10) $\begin{array}{rl}& \beta =0.4,{\alpha }_1=0.5,{\alpha }_2=0.1,\rho =0.4,\sigma =0.25,\\ & \kappa =0.5,{\gamma }_1=0.14,{\gamma }_2=0.14,{\gamma }_3=0.25,\mu =2.3\times {10}^{-4}.\end{array}$(10) The corresponding basic reproduction number is ${\mathcal{R}}_0={\mathcal{R}}_1+{\mathcal{R}}_2+{\mathcal{R}}_3\approx 1.88$ with ${\mathcal{R}}_1=0.12,{\mathcal{R}}_2=1.71$ and ${\mathcal{R}}_3=0.05$. So, the percentages of contribution of the infected states $I_s,I_a$ and Q to ${\mathcal{R}}_0$ are 6.62%, 90.73% and 2.65%, respectively. This suggests that asymptomatic transmission is a leading cause of HFMD persistence. We then calculate their sensitivity indices as shown in Figure 5(a). Obviously, ${\mathcal{R}}_0$ is most sensitive to the transmission coefficient, β, followed by the recovery rate of asymptomatic infections, ${\gamma }_2$, and the proportion of infections being symptomatic, ρ, and is least sensitive to the progression rate from the exposed state to the infectious state, σ, the progression rate leaving children group, μ, and the recovery rate of the symptomatic people, ${\gamma }_1$.

Figure 5. Sensitivity analysis for ${\mathcal{R}}_0$ with respect to model parameters: (a) sensitivity indices (SI), and (b) partial rank correlation coefficients (PRCC).

The conclusion from the sensitivity indices can be biased since it only measures the influence of a single parameter on model output which could be strongly relied on the choice of a specific parameter set. To improve the robustness of sensitivity analysis, we investigate the impact of global parameter variations on model outcome. By using the Latin Hypercube Sampling (LHS) method [29], we generate ${10}^6$ random parameter sets with ranges in Table 1. Each input parameter is assumed to be uniformly distributed. Then we calculate the partial rank correlation coefficient (PRCC) of ${\mathcal{R}}_0$ with respect to each involved parameter (see Figure 5(b)). The global sensitivity analysis gives a similar conclusion, whereas the main difference is that the relative contact rate of the quarantined state, ${\alpha }_2$, is moderately positively correlated to ${\mathcal{R}}_0$. In reality, it is hard to detect and treat asymptomatic cases and change the symptomatic ratio. This indicates that reducing contacts between susceptible and infectious individuals and implementing strict quarantine are vital to HFMD control.

4.3. Numerical simulations

Example 4.1

Behaviour change on infection size

We notice that the behaviour change function $f(I_s,Q)$ does not influence ${\mathcal{R}}_0$. Thus, variable behaviour change has no impact on disease persistence. However, it follows from Remark 3.4 that the change can protect some susceptibles from contracting the disease. To quantitatively examine the role of behaviour change on a HFMD epidemic wave, we consider $\mathrm{\Lambda }=463$ and the remaining parameter setting is the same as (10(10) $\begin{array}{rl}& \beta =0.4,{\alpha }_1=0.5,{\alpha }_2=0.1,\rho =0.4,\sigma =0.25,\\ & \kappa =0.5,{\gamma }_1=0.14,{\gamma }_2=0.14,{\gamma }_3=0.25,\mu =2.3\times {10}^{-4}.\end{array}$(10) ). The initial conditions are as follows: $(S(0),E(0),I_s(0),I_a(0),Q(0),R(0))=(1.11\times {10}^6,10,0,0,0,9\times {10}^5).$Again the basic reproduction number is ${\mathcal{R}}_0\approx 1.88>1$. We adopt the same behaviour change function as given in (9(9) $f(I_s,Q)={\displaystyle \frac{1}{1+m_1{I}_s+m_{2}Q}},$(9) ) and select three scenarios with different level of behaviour change (BC):

1. no behaviour change: $m_1=m_2=0$;

2. weak behaviour change: $m_1=6\times {10}^{-4}$ and $m_2=3\times {10}^{-4}$;

3. strong behaviour change: $m_1=9\times {10}^{-3}$ and $m_2=4.5\times {10}^{-3}$.

Figure 6(a) shows the total number of symptomatically infected individuals, i.e. $I_s+Q$, with $m_1$ and $m_2$ increasing from zero to baseline values to 15-fold of them. We can see that the introduction of infectives always leads to an outbreak and the disease spread attains its first peak almost at the same time under no and weak behaviour change. The endemic equilibria ${\mathit{E}}^{\ast }=(S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast },R^{\ast })$ and the nonsusceptible ratios $1-\frac{S^{\ast }}{N^{\ast }}$ associated to the three scenarios are: $\begin{array}{rl}& (1.068\times {10}^6,869,136,929,271,9.427\times {10}^5)\text{and}46.9\mathrm{\%},\text{no behaviour change,}\\ & (1.202\times {10}^6,746,116,798,233,8.094\times {10}^5)\text{and}40.3\mathrm{\%},\text{weak behaviour change,}\\ & (1.731\times {10}^6,260,41,278,81,2.819\times {10}^5)\text{and}14.0\mathrm{\%},\text{strong behaviour change,}\end{array}$respectively. Weak behaviour change can sharply reduce the peak size (the maximum number of symptomatic infections) but moderately lower the endemic level, while strong behaviour change reduces both peak size and endemic level dramatically. When there is no behaviour change, the solution converges to the endemic equilibrium in an oscillating way. Stronger behaviour change results in less or no damped oscillations over time. Figure 6(b) illustrates that the healthy people eventually reduce their contacts by 12% and 42% through weak and strong behaviour change, respectively. This can also be obtained by noting $f(I_s(t),Q(t))\to f(I_s^{\ast },Q^{\ast })$ as $t\to \mathrm{\infty }$. In particular, under strong behaviour change, 42% reduction in the contact rate of healthy people can lead to 70% reduction in symptomatic infections. The respective percentages of contribution of $I_s,I_a$ and Q to new infections change from 6.62%, 90.73% and 2.65% (same as ${\mathcal{R}}_0$) to 7.45%, 89.57% and 2.98% to 10.73%, 84.98% and 4.29% as behaviour change scenario changes from (a) to (b) to (c).

Figure 6. Simulations of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q,R)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q,R)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& ={\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q-\mu R,\end{array}$(1) ) for (a) the square root of the total number of symptomatic infections, $\sqrt{I_s+Q}$, and (b) the behaviour change function, $f(I_s,Q)$, versus time under no behaviour change (red dashed line), weak behaviour change (black solid line), and strong behaviour change (blue dotted line). See main text for parameter values and initial condition.

Furthermore, using the same parameter setting and initial condition, the dependences of the total symptomatic infections and the relative contact rate of the healthy people at the endemic equilibrium, i.e. $I_s^{\ast }+Q^{\ast }$ and $f(I_s^{\ast },Q^{\ast })$, against the behaviour change parameters, $m_1$ and $m_2$, are plotted in Figure 7. Clearly, the more sensitive the public is to an increase in the number of new or recent infections, the lower the number of cases and contacts. Note that the contours of $I_s^{\ast }+Q^{\ast }$ and $f(I_s^{\ast },Q^{\ast })$ are straight lines in the $m_1$-$m_2$ plane. Indeed, given behaviour change parameters ${\tilde{m}}_1\ge 0$ and ${\tilde{m}}_2\ge 0$, denote the corresponding endemic equilibrium by ${\tilde{\mathit{E}}}^{\ast }=({\tilde{S}}^{\ast },{\tilde{E}}^{\ast },{\tilde{I}}_s^{\ast },{\tilde{I}}_a^{\ast },{\tilde{Q}}^{\ast },{\tilde{R}}^{\ast })$, then $f({\tilde{I}}_s^{\ast },{\tilde{Q}}^{\ast })=\frac{1}{1+{\tilde{m}}_1{\tilde{I}}_s^{\ast }+{\tilde{m}}_2{\tilde{Q}}^{\ast }}\text{and}\kappa {\tilde{I}}_s^{\ast }=({\gamma }_3+\mu ){\tilde{Q}}^{\ast }.$Therefore, for any pair of parameter values on the line segment $\mathcal{L}:\{(m_1,m_2)\in {\mathbb{R}}_{+}^2:m_1{\tilde{I}}_s^{\ast }+m_2{\tilde{Q}}^{\ast }={\tilde{m}}_1{\tilde{I}}_s^{\ast }+{\tilde{m}}_2{\tilde{Q}}^{\ast }\},$the associated unique endemic equilibrium remains unchanged. More specifically, the line is $(m_1-{\tilde{m}}_1){\tilde{I}}_s^{\ast }+(m_2-{\tilde{m}}_2){\tilde{Q}}^{\ast }=0,\text{i.e.}({\gamma }_3+\mu )(m_1-{\tilde{m}}_1)+\kappa (m_2-{\tilde{m}}_2)=0.$Hence, the two quantities $I_s^{\ast }+Q^{\ast }$ and $f(I_s^{\ast },Q^{\ast })$ are constant on $\mathcal{L}$.

Example 4.2

Quarantine on infection size

We now examine the impact of quarantine on symptomatic cases. Note that quarantine strategy is characterized by κ and ${\alpha }_2$, which are related to the speed and efficiency of quarantine. There is yet no specific treatment for HFMD, so we assume that quarantine does not accelerate the recovery process.

Figure 7. The contour plots of the total number of symptomatic infections and the behaviour change function at the endemic equilibrium, $I_s^{\ast }+Q^{\ast }$ and $f(I_s^{\ast },Q^{\ast })$, versus parameters $m_1$ and $m_2$. See main text for parameter values.

Consider the same set of parameter values as given in Example 4.1 with the exception of κ and ${\alpha }_2$. The initial conditions of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q,R)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q,R)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& ={\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q-\mu R,\end{array}$(1) ) are $(S(0),E(0),I_s(0),I_a(0),Q(0),R(0))=(2.01\times {10}^6,10,0,0,0,0).$For fixed ${\alpha }_2=0.1$, the total numbers of symptomatic cases $I_s+Q$ over time are plotted in Figure 8(a) under the three behaviour change scenarios (a)–(c) as κ increases from 0.1 to 0.3 to 0.9. The basic reproduction numbers ${\mathcal{R}}_0$ correspond to $\kappa =0.1,0.3$ and 0.9 are 2.07, 1.93 and 1.84, respectively. When there is no change in behaviour, fast quarantine reduces the peak size and the epidemic size of symptomatic infections considerably and delays the peak time (the time of the maximum number of infected humans) by 2-4 weeks. In the presence of behaviour change, fast quarantine only has a mild effect on disease propagation. Interestingly, it suggests that a higher but shorter epidemic (red curves) may have a larger epidemic size than a lower but much longer epidemic (blue or black curves). Meanwhile, given $\kappa =0.3$, Figure 8(b) shows how the total numbers of symptomatic infections evolve under different levels of behaviour change as the value of ${\alpha }_2$ decreases from 0.4 to 0.1 to 0.025. So, improving the efficiency of quarantine plays a similar but weaker role in containing disease spread which means that timely quarantine is more important than strict quarantine. It is worth mentioning that quarantine significantly affects the infection size even though it has little impact on ${\mathcal{R}}_0$.

Figure 8. Simulations of model (1(1) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q,R)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q,R)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\\ \frac{\mathrm{d}R}{\mathrm{d}t}& ={\gamma }_1{I}_s+{\gamma }_2{I}_a+{\gamma }_{3}Q-\mu R,\end{array}$(1) ) for the cube root of the total number of symptomatic infections, $\sqrt[3]{I_s+Q}$, under different behaviour change and (a) quarantine rate $\kappa =0.1$, 0.3 and 0.9, and (b) relative contact rate for the quarantined individuals ${\alpha }_2=0.4$, 0.1 and 0.025.

5. Discussion

HFMD is a childhood contagious disease that spreads widely in China. In recent years, a number of mathematical modelling studies on HFMD have been done with the consideration of seasonal variation in contact pattern, vaccine introduction, environmental contamination and so on [23, 34, 43]. During a disease outbreak, changes in behaviour driven by governmental action or individual reaction are very common. Their role in disease control and elimination has been increasingly considered by modellers [12, 19, 30]. In this paper, we proposed a deterministic endemic model with imperfect quarantine to explore the impact of behaviour change on the spread of HFMD. The total population is split into six compartments where the symptomatic and quarantined compartments obey constant behaviour change and the behaviour change of other compartments depends on the number of newly reported cases and recently cumulative cases. We defined a biologically feasible region and derived the basic reproduction number ${\mathcal{R}}_0$ of the model. Using the Lyapunov functional method and the persistence theory, a threshold dynamic result was obtained, i.e. the disease-free equilibrium is globally asymptotically stable as ${\mathcal{R}}_0\le 1$ and the disease is uniformly persistent otherwise. Moreover, there exists a unique endemic equilibrium if ${\mathcal{R}}_0>1$ and a sufficient condition for the global stability of the endemic equilibrium was given.

In addition, we fitted our model to the reported HFMD case data in Shanghai and found that the basic reproduction number is around 2.41. Based on the fitting result and parameter ranges, we conducted both local and global sensitivity analysis for ${\mathcal{R}}_0$ in terms of model parameters. Control measures that reduce close contacts and improve quarantine adherence can be successful in constraining HFMD. Two numerical examples were presented to analyse the influence of behaviour change and quarantine on reducing disease burden. In the first example, we compared the total numbers of symptomatic infections under different behaviour changes. Although incidence-based or prevalence-based behaviour change of healthy people does not affect the disease persistence, it substantially lowers the endemic level and flattens the curve so that the health care capacity can meet the demands during an outbreak. Ignoring behaviour change can potentially overestimate the disease magnitude and the contribution of asymptomatic transmission. The stronger the behaviour change to an epidemic, the smaller the epidemic size. Thus, it is beneficial to run more health education campaigns during high transmission season. In the second example, we demonstrated that fast quarantine of symptomatic cases can significantly curtail the infection size and delay the peak time in the absence of behaviour change. Nevertheless, behaviour change may weaken the benefit of quarantine.

To some extent, we generalize the HFMD model of Ma et al. [27] by considering behaviour change and imperfect quarantine. Although our modelling and analysis are based on HFMD, they are applicable to some other infectious diseases such as measles and mumps. Most existing epidemic models with behaviour change either assume that all individuals have the same behaviour change depending on the number of infected individuals [7, 10, 42], or divide the total population into multiple groups and different groups have different but constant behaviour change [5, 30, 46]. The main novelty of our model is that some states/individuals have variable behaviour change while others have constant behaviour change. Meanwhile, we assume that the variable behaviour change is determined by the number of confirmed cases and its change rate. The idea is somewhat similar to that of Xiao et al. [45] but we use the number of symptomatic cases to replace the change rate of infected people which facilitates the mathematical analysis. A related form of contact rate was introduced by Liu et al. [24] but the epidemiological meaning is different. We found that the transmission from mild or asymptomatic patients plays a dominant role in maintaining HFMD spreading. This result is surprising and arguable [34]. On the one hand, symptomatic cases of HFMD are often easy to identify, sharply reducing their transmission potential. On the other hand, the presumed high asymptomatic ratio and the equal infectivity for the asymptomatic and symptomatic individuals may exaggerate the contribution of asymptomatic transmission.

There are many possibilities to improve and generalize the current work. The global asymptotic stability of the endemic equilibrium is generally unclear. Fitting the model to long term rather than one-year data is challenging but more convincing. Seasonal epidemics of HFMD indicate that a periodic epidemic model with time-dependent transmission coefficient may be favourable [23]. Some additional biological factors like environment-to-human transmission route [42], coexistence of multiple enteroviruses [32], EV71 vaccination [34], population movement [16], waning immunity, age structure, stochasticity, and varying population size can be added to the model. A full understanding of how the emergency response plan is implemented in schools and kindergartens can help us to build a more realistic model. How to mechanistically determine and quantify the behaviour change function? To what extent can behaviour change reduce the final size of an epidemic as vital dynamics are ignored (see Brauer [6])? Comparing the constant behaviour change, variable behaviour change and mixed behaviour change approaches with epidemiological and behavioural data is necessary for the selection and application of behaviour change models. It is desirable to develop a general modelling framework with heterogeneities in behaviour change. Since ${\mathcal{R}}_0$ is defined at the disease-free equilibrium at which there is no behaviour change, the reproduction number or disease dynamics are unaffected by behaviour change. In reality, long-term behaviour change may eradicate an infectious disease in a specific region. One way to construct a behaviour change function that still works when the initial conditions are near the disease-free space is to assume that the contact rate depends on the cumulative number of cases [17]. Changes in behaviour are not limited to contact rate but also include vaccination rate, recruitment rate [40], movement rate [26] and so on. Besides the case and death counts, some other factors can strengthen or weaken behaviour change. The ongoing COVID-19 pandemic indicates that individual response to epidemic waves varies significantly with vaccine development and distribution. The biennial cycle of HFMD incidence indicates that there should be more cases in 2020 than in 2019. However, due to the COVID-19 outbreak, there were 5,585 reported cases of HFMD in 2020 compared to 24,615 cases in 2019 in Shanghai while there were 761,355 reported HFMD cases in 2020 compared to 1,918,830 cases in 2019 in China [1, 49]. It is attractive to understand the transmission dynamics of HFMD under the COVID-19-induced behaviour change. We must realize that behaviour change like limiting contact and travel could be costly and unsustainable. A cost-effectiveness analysis of behaviour change measures is indispensable. In the era of mobile internet and artificial intelligence, the rapid development in data collection, storage, and extraction methods makes behaviour change models for infectious diseases more applicable but there are still quite a few challenges to overcome [13].

Acknowledgments

We thank the two referees for their valuable comments, Dr. Jim Cushing for handling the submission, and Drs. Jifa Jiang and Yijun Lou for helpful suggestions on the draft.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Natural Science Foundation of Shanghai [grant numbers 20ZR1440600 and 20JC1413800], the National Natural Science Foundation of China [grant number 12071300], the Key Research Projects of Beijing Natural Science Foundation-Haidian District Joint Fund [grant number L192012], and the CSU Office of Research through a startup grant.

Appendix. Proof of Theorem 3.5

Proof.

Let ${\mathit{E}}^{\ast }=(S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })$ denote the endemic equilibrium of model (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) which satisfies (5). Define $V_1=S^{\ast }g\left(\frac{S}{S^{\ast }}\right)+E^{\ast }g\left(\frac{E}{E^{\ast }}\right),V_2=I_s^{\ast }g\left(\frac{I_s}{I_s^{\ast }}\right),V_3=I_a^{\ast }g\left(\frac{I_a}{I_a^{\ast }}\right),V_4=Q^{\ast }g\left(\frac{Q}{Q^{\ast }}\right),$where $g(x)=x-1-\ln x\ge 0$ with equality if and only if x = 1.

Consider a Lyapunov function $V=l_1{V}_1+l_2{V}_2+l_3{V}_3+l_4{V}_4,$where positive constants $l_i$, i = 1, 2, 3, 4 are to be determined. The derivative of V along the solution of system (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) is $V^{\prime }=l_1(1-{\displaystyle \frac{S^{\ast }}{S}})S^{\prime }+l_1(1-{\displaystyle \frac{E^{\ast }}{E}})E^{\prime }+l_2(1-{\displaystyle \frac{I_s^{\ast }}{I_s}})I_s^{\prime }+l_3(1-{\displaystyle \frac{I_a^{\ast }}{I_a}})I_a^{\prime }+l_4(1-{\displaystyle \frac{Q^{\ast }}{Q}})Q^{\prime }.$For convenience, we write $f=f(I_s,Q)$ and $f^{\ast }=f(I_s^{\ast },Q^{\ast })$. Using (5a(5a) $\begin{array}{rl}& \mathrm{\Lambda }-\lambda (S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })S^{\ast }-\mu S^{\ast }=0,\end{array}$(5a) ) and (5b(5b) $\begin{array}{rl}& \lambda (S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })S^{\ast }-(\sigma +\mu )E^{\ast }=0,\end{array}$(5b) ), $\begin{array}{rl}{V}_1^{\prime }& =(1-{\displaystyle \frac{S^{\ast }}{S}})(\beta {\displaystyle \frac{{\alpha }_1{I}_s^{\ast }+f^{\ast }{I}_a^{\ast }+{\alpha }_2{Q}^{\ast }}{N^{\ast }}}{S}^{\ast }+\mu S^{\ast }-\beta {\displaystyle \frac{{\alpha }_1{I}_s+f{I}_a+{\alpha }_{2}Q}{N^{\ast }}}S-\mu S)\\ & +(1-{\displaystyle \frac{E^{\ast }}{E}})(\beta {\displaystyle \frac{{\alpha }_1{I}_s+f{I}_a+{\alpha }_{2}Q}{N^{\ast }}}S-\beta {\displaystyle \frac{{\alpha }_1{I}_s^{\ast }+f^{\ast }{I}_a^{\ast }+{\alpha }_2{Q}^{\ast }}{N^{\ast }}}{S}^{\ast }\cdot {\displaystyle \frac{E}{E^{\ast }}})\\ & =-{\displaystyle \frac{\mu (S-S^{\ast }{)}^2}{S}}+\frac{{\alpha }_1\beta I_s^{\ast }{S}^{\ast }}{N^{\ast }}(2+{\displaystyle \frac{I_s}{I_s^{\ast }}}-{\displaystyle \frac{S^{\ast }}{S}}-{\displaystyle \frac{E}{E^{\ast }}}-{\displaystyle \frac{S{I}_s{E}^{\ast }}{S^{\ast }{I}_s^{\ast }E}})\\ & +\frac{f^{\ast }\beta I_a^{\ast }{S}^{\ast }}{N^{\ast }}(2+{\displaystyle \frac{f{I}_a}{f^{\ast }{I}_a^{\ast }}}-{\displaystyle \frac{S^{\ast }}{S}}-{\displaystyle \frac{E}{E^{\ast }}}-{\displaystyle \frac{fS{I}_a{E}^{\ast }}{f^{\ast }{S}^{\ast }{I}_a^{\ast }E}})\\ & +\frac{{\alpha }_2\beta Q^{\ast }{S}^{\ast }}{N^{\ast }}(2+{\displaystyle \frac{Q}{Q^{\ast }}}-{\displaystyle \frac{S^{\ast }}{S}}-{\displaystyle \frac{E}{E^{\ast }}}-{\displaystyle \frac{SQ{E}^{\ast }}{S^{\ast }{Q}^{\ast }E}}).\end{array}$Then we have $\begin{array}{rl}& 2+{\displaystyle \frac{I_s}{I_s^{\ast }}}-{\displaystyle \frac{S^{\ast }}{S}}-{\displaystyle \frac{E}{E^{\ast }}}-{\displaystyle \frac{S{I}_s{E}^{\ast }}{S^{\ast }{I}_s^{\ast }E}}=(1-{\displaystyle \frac{S{I}_s{E}^{\ast }}{S^{\ast }{I}_s^{\ast }E}})+(1-{\displaystyle \frac{S^{\ast }}{S}})+{\displaystyle \frac{I_s}{I_s^{\ast }}}-{\displaystyle \frac{E}{E^{\ast }}}\\ & \le \ln {\displaystyle \frac{S^{\ast }{I}_s^{\ast }E}{S{I}_s{E}^{\ast }}}+\ln \frac{S}{S^{\ast }}+{\displaystyle \frac{I_s}{I_s^{\ast }}}-{\displaystyle \frac{E}{E^{\ast }}}=\ln \frac{I_s^{\ast }E}{I_s{E}^{\ast }}+{\displaystyle \frac{I_s}{I_s^{\ast }}}-{\displaystyle \frac{E}{E^{\ast }}}\\ & =\ln {\displaystyle \frac{E}{E^{\ast }}}-\ln {\displaystyle \frac{I_s}{I_s^{\ast }}}+{\displaystyle \frac{I_s}{I_s^{\ast }}}-{\displaystyle \frac{E}{E^{\ast }}},\end{array}$where the equality holds if and only if $S^{\ast }=S$ and $\frac{E}{E^{\ast }}=\frac{I_s}{I_s^{\ast }}$. Similarly, $2+{\displaystyle \frac{Q}{Q^{\ast }}}-{\displaystyle \frac{S^{\ast }}{S}}-{\displaystyle \frac{E}{E^{\ast }}}-{\displaystyle \frac{SQ{E}^{\ast }}{S^{\ast }{Q}^{\ast }E}}\le \ln {\displaystyle \frac{E}{E^{\ast }}}-\ln {\displaystyle \frac{Q}{Q^{\ast }}}+{\displaystyle \frac{Q}{Q^{\ast }}}-{\displaystyle \frac{E}{E^{\ast }}}$with equality if and only if $S^{\ast }=S$ and $\frac{E}{E^{\ast }}=\frac{Q}{Q^{\ast }}$. By assumption, $\begin{array}{rl}& 2+{\displaystyle \frac{f{I}_a}{f^{\ast }{I}_a^{\ast }}}-{\displaystyle \frac{S^{\ast }}{S}}-{\displaystyle \frac{E}{E^{\ast }}}-{\displaystyle \frac{fS{I}_a{E}^{\ast }}{f^{\ast }{S}^{\ast }{I}_a^{\ast }E}}\\ & =-(1-\frac{f{I}_a}{f^{\ast }{I}_a^{\ast }})(1-\frac{f^{\ast }}{f})+3-\frac{S^{\ast }}{S}-{\displaystyle \frac{fS{I}_a{E}^{\ast }}{f^{\ast }{S}^{\ast }{I}_a^{\ast }E}}-\frac{f^{\ast }}{f}-\frac{E}{E^{\ast }}+\frac{I_a}{I_a^{\ast }}\\ & \le -(\frac{S^{\ast }}{S}-1)-({\displaystyle \frac{fS{I}_a{E}^{\ast }}{f^{\ast }{S}^{\ast }{I}_a^{\ast }E}}-1)-(\frac{f^{\ast }}{f}-1)-\frac{E}{E^{\ast }}+\frac{I_a}{I_a^{\ast }}\\ & \le -\ln (\frac{S^{\ast }}{S}\cdot {\displaystyle \frac{fS{I}_a{E}^{\ast }}{f^{\ast }{S}^{\ast }{I}_a^{\ast }E}}\cdot \frac{f^{\ast }}{f})-\frac{E}{E^{\ast }}+\frac{I_a}{I_a^{\ast }}\\ & =\ln {\displaystyle \frac{E}{E^{\ast }}}-\ln {\displaystyle \frac{I_a}{I_a^{\ast }}}+{\displaystyle \frac{I_a}{I_a^{\ast }}}-{\displaystyle \frac{E}{E^{\ast }}}\end{array}$with equality if and only if $S^{\ast }=S$, $\frac{E}{E^{\ast }}=\frac{I_a}{I_a^{\ast }}$ and $f(I_s,Q)=f(I_s^{\ast },Q^{\ast })$. Therefore, $\begin{array}{rl}{V}_1^{\prime }& \le \frac{{\alpha }_1\beta I_s^{\ast }{S}^{\ast }}{N^{\ast }}(\ln {\displaystyle \frac{E}{E^{\ast }}}-\ln {\displaystyle \frac{I_s}{I_s^{\ast }}}+{\displaystyle \frac{I_s}{I_s^{\ast }}}-{\displaystyle \frac{E}{E^{\ast }}})+\frac{f^{\ast }\beta I_a^{\ast }{S}^{\ast }}{N^{\ast }}(\ln {\displaystyle \frac{E}{E^{\ast }}}-\ln {\displaystyle \frac{I_a}{I_a^{\ast }}}+{\displaystyle \frac{I_a}{I_a^{\ast }}}-{\displaystyle \frac{E}{E^{\ast }}})\\ & +\frac{{\alpha }_2\beta Q^{\ast }{S}^{\ast }}{N^{\ast }}(\ln {\displaystyle \frac{E}{E^{\ast }}}-\ln {\displaystyle \frac{Q}{Q^{\ast }}}+{\displaystyle \frac{Q}{Q^{\ast }}}-{\displaystyle \frac{E}{E^{\ast }}}).\end{array}$Similarly, by using (5c(5c) $\begin{array}{rl}& \rho \sigma E^{\ast }-(\kappa +{\gamma }_1+\mu )I_s^{\ast }=0,\end{array}$(5c) ), $\begin{array}{rl}{V}_2^{\prime }& =(1-{\displaystyle \frac{I_s^{\ast }}{I_s}})(\rho \sigma E-{\displaystyle \frac{\rho \sigma E^{\ast }}{I_s^{\ast }}}{I}_s)=\rho \sigma E^{\ast }(1+{\displaystyle \frac{E}{E^{\ast }}}-{\displaystyle \frac{I_s}{I_s^{\ast }}}-{\displaystyle \frac{I_s^{\ast }E}{I_s{E}^{\ast }}})\\ & \le \rho \sigma E^{\ast }({\displaystyle \frac{E}{E^{\ast }}}-{\displaystyle \frac{I_s}{I_s^{\ast }}}-\ln {\displaystyle \frac{E}{E^{\ast }}}+\ln {\displaystyle \frac{I_s}{I_s^{\ast }}})\end{array}$with equality if and only if $\frac{E}{E^{\ast }}=\frac{I_s}{I_s^{\ast }}$. By using (5d(5d) $\begin{array}{rl}& (1-\rho )\sigma E^{\ast }-({\gamma }_2+\mu )I_a^{\ast }=0,\end{array}$(5d) ), $\begin{array}{rl}{V}_3^{\prime }& =(1-{\displaystyle \frac{I_a^{\ast }}{I_a}})((1-\rho )\sigma E-{\displaystyle \frac{(1-\rho )\sigma E^{\ast }}{I_a^{\ast }}}{I}_a)=(1-\rho )\sigma E^{\ast }(1+{\displaystyle \frac{E}{E^{\ast }}}-{\displaystyle \frac{I_a}{I_a^{\ast }}}-{\displaystyle \frac{I_a^{\ast }E}{I_a{E}^{\ast }}})\\ & \le (1-\rho )\sigma E^{\ast }({\displaystyle \frac{E}{E^{\ast }}}-{\displaystyle \frac{I_a}{I_a^{\ast }}}-\ln {\displaystyle \frac{E}{E^{\ast }}}+\ln {\displaystyle \frac{I_a}{I_a^{\ast }}})\end{array}$with equality if and only if $\frac{E}{E^{\ast }}=\frac{I_a}{I_a^{\ast }}$. By using (5e(5e) $\begin{array}{rl}& \kappa I_s^{\ast }-({\gamma }_3+\mu )Q^{\ast }=0.\end{array}$(5e) ), $\begin{array}{rl}{V}_4^{\prime }& =(1-{\displaystyle \frac{Q^{\ast }}{Q}})(\kappa I_s-{\displaystyle \frac{\kappa I_s^{\ast }}{Q^{\ast }}}Q)=\kappa I_s^{\ast }(1+{\displaystyle \frac{I_s}{I_s^{\ast }}}-{\displaystyle \frac{Q}{Q^{\ast }}}-{\displaystyle \frac{Q^{\ast }{I}_s}{Q{I}_s^{\ast }}})\\ & \le \kappa I_s^{\ast }({\displaystyle \frac{I_s}{I_s^{\ast }}}-{\displaystyle \frac{Q}{Q^{\ast }}}-\ln {\displaystyle \frac{I_s}{I_s^{\ast }}}+\ln {\displaystyle \frac{Q}{Q^{\ast }}})\end{array}$with equality if and only if $\frac{I_s}{I_s^{\ast }}=\frac{Q}{Q^{\ast }}$.

Combining the above estimates gives $\begin{array}{rl}{V}^{\prime }& \le (-{\displaystyle \frac{l_1\beta S^{\ast }}{N^{\ast }}}({\alpha }_1{I}_s^{\ast }+f^{\ast }{I}_a^{\ast }+{\alpha }_2{Q}^{\ast })+l_2\rho \sigma E^{\ast }+l_3(1-\rho )\sigma E^{\ast })({\displaystyle \frac{E}{E^{\ast }}}-\ln {\displaystyle \frac{E}{E^{\ast }}})\\ & +(l_1\frac{{\alpha }_1\beta I_s^{\ast }{S}^{\ast }}{N^{\ast }}-l_2\rho \sigma E^{\ast }+l_4\kappa I_s^{\ast })({\displaystyle \frac{I_s}{I_s^{\ast }}}-\ln {\displaystyle \frac{I_s}{I_s^{\ast }}})\\ & +(l_1\frac{f^{\ast }\beta I_a^{\ast }{S}^{\ast }}{N^{\ast }}-l_3(1-\rho )\sigma E^{\ast })({\displaystyle \frac{I_a}{I_a^{\ast }}}-\ln {\displaystyle \frac{I_a}{I_a^{\ast }}})\\ & +(l_1\frac{{\alpha }_2\beta Q^{\ast }{S}^{\ast }}{N^{\ast }}-l_4\kappa I_s^{\ast })({\displaystyle \frac{Q}{Q^{\ast }}}-\ln {\displaystyle \frac{Q}{Q^{\ast }}}).\end{array}$Setting the coefficient of each term on the right-hand side of the inequality to zero leads to a system of linear equations with respect to $l_1,l_2,l_3$ and $l_4$, from which we obtain $\begin{array}{rl}{l}_1& =\rho (1-\rho )\sigma E^{\ast },l_2={\displaystyle \frac{(1-\rho )\beta S^{\ast }({\alpha }_1{I}_s^{\ast }+{\alpha }_2{Q}^{\ast })}{N}},\\ l_3& ={\displaystyle \frac{\rho \beta f(I_s^{\ast },Q^{\ast })S^{\ast }{I}_a^{\ast }}{N}},l_4={\displaystyle \frac{\rho (1-\rho )\sigma {\alpha }_2\beta S^{\ast }{E}^{\ast }{Q}^{\ast }}{\kappa N{I}_s^{\ast }}}.\end{array}$Under this choice of $l_i$ for i = 1, 2, 3, 4, we have $V^{\prime }\le 0$. Furthermore, $V^{\prime }=0$ implies that $S=S^{\ast },f(I_s,Q)=f(I_s^{\ast },Q^{\ast }),E=q{E}^{\ast },I_s=q{I}_s^{\ast },I_a=q{I}_a^{\ast }\text{and}Q=q{Q}^{\ast }$for some q>0. Substituting these relations into the first equation of (2(2) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =\mathrm{\Lambda }-\lambda (S,E,I_s,I_a,Q)S-\mu S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}& =\lambda (S,E,I_s,I_a,Q)S-(\sigma +\mu )E,\\ \frac{\mathrm{d}{I}_s}{\mathrm{d}t}& =\rho \sigma E-(\kappa +{\gamma }_1+\mu )I_s,\\ \frac{\mathrm{d}{I}_a}{\mathrm{d}t}& =(1-\rho )\sigma E-({\gamma }_2+\mu )I_a,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}& =\kappa I_s-({\gamma }_3+\mu )Q,\end{array}$(2) ) gives $\begin{array}{rl}0& =\mathrm{\Lambda }-\lambda (S^{\ast },q{E}^{\ast },q{I}_s^{\ast },q{I}_a^{\ast },q{Q}^{\ast })S^{\ast }-\mu S^{\ast }=\mathrm{\Lambda }-\lambda (S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })S^{\ast }-\mu S^{\ast }\\ & \iff \lambda (S^{\ast },q{E}^{\ast },q{I}_s^{\ast },q{I}_a^{\ast },q{Q}^{\ast })=\lambda (S^{\ast },E^{\ast },I_s^{\ast },I_a^{\ast },Q^{\ast })\\ & \iff \beta {\displaystyle \frac{{\alpha }_{1}q{I}_s^{\ast }+f^{\ast }q{I}_a^{\ast }+{\alpha }_{2}q{Q}^{\ast }}{N^{\ast }}}=\beta {\displaystyle \frac{{\alpha }_1{I}_s^{\ast }+f^{\ast }{I}_a^{\ast }+{\alpha }_2{Q}^{\ast }}{N^{\ast }}}\iff q=1.\end{array}$Thus, $V^{\prime }=0$ holds only at ${\mathit{E}}^{\ast }$. By the standard Lyapunov stability theorem, the global stability of ${\mathit{E}}^{\ast }$ is established.