A multi-group model for estimating the transmission rate of hand, foot and mouth disease in mainland China

1.   Introduction

Hand, foot and mouth disease (HFMD) is an emerging illness which usually affects infants and children by coxsackievirus A16 (CAV 16) and human enterovirus 71 (EV71). Besides, many other strains of coxsackievirus and enterovirus are responsible for spreading HFMD as well. The majority of infected individuals are children under 5 years and the ones most likely to develop this disease are young children under 3 years, but this illness may be transmitted among adults ^[1]. Because of less immune and self-aware, children are more susceptible to infection than adults. Typical symptom of HFMD is fever (37.8 ℃ to 38.9 ℃), painful sores in the mouth and a rash with blisters on the hands, feet and buttocks. Diagnosis is usually made only through signs and symptoms. If the diagnosis is not clear, a throat swab or stool specimen can be taken to identify the virus by culture ^[2]. Generally, the incubation period is 2-7 days. HFMD is highly contagious and can be transmitted through nasopharyngeal secretions such as saliva or nasal mucus, direct contact or fecal-oral transmission ^[3]. Even though few infected children and most infected adults have no symptoms, they are contagious since they can transmit virus, who are subclinical cases.

HFMD is most common in mainland China and tends to break out in spring and autumn especially ^[1]. Preventive measures include avoiding direct contact with infected persons (such as keep the infected children away from childcare or school), proper cleaning of shared equipment, disinfection of contaminated surfaces and appropriate hand hygiene. These measures have proven to be effective in reducing viral transmission responsible for HFMD ^[4,5]. Many affected countries have adopted routine control measures similar to pandemic preparedness plans, including surveillance, mandatory reporting, isolation, school closure and social alienation. A vaccine known as the EV71 vaccine is available to prevent HFMD in China as of December 2015 ^[6,8], but the vaccine does not belong to National Immunization Program by Ministry of the People's Republic of China, perhaps because of safety concerns, and is also not widely used in hospital. At present, people only in a few countries including China are vaccinated against HFMD ^[9]. So for now, it is incredible that there is not specific curative treatment for HFMD ^[10].

HFMD usually doesn't require medication and it can resolve itself. Disease control usually focuses on alleviating symptoms, and pain from the sores may be eased with the use of analgesic medications. In most cases, the disease is mild and self-limiting ^[11], but the more serious clinical symptoms are neurological abnormalities and even other serious complications, such as myocarditis, pulmonary edema and aseptic meningoencephalitis in few children. Some severely affected patients may die due to agammaessive malignancy of the disease ^[2].

Epidemiological models have become important tools in understanding well the spread and control of infectious diseases ^[7]. Recently, there are several types of mathematical models that have been used to investigate the transmission dynamics and predict HFMD infections ^{[8,28,34,35,36,37,38,41]}, and the estimates of basic reproduction number can be seen in Table 1. It may be meaningful to consider a model emerged HFMD cases at different geographical locations, ages or other categories, then identify the high-risk group and main population of transmission. Moreover, taking different internal structures of the host population and the transmission properties of infectious diseases into account, a heterogeneous host population can be partitioned into several homogeneous subetaoups, according to various characteristics of individuals, such as age, contact patterns, social and economic status, profession and demographical distribution ^[12,13,14]. This is known as a multi-group model. One of the earliest multi-group models was proposed by Lajmanovich and Yorke ^[12] for the transmission of gonorrhea.

Obviously, in the multi-group model, we have to consider the interactions within a subetaoup as well as among different subetaoups in the course of the transmission of infectious diseases. It can produce more interesting and complicated scenarios of disease transmission in the multi-group model. A tremendous variety of multi-group models have been formulated, analyzed, and applied to many infectious diseases, see ^{[15,16,17,18,19,20,23]} for example. Multi-group epidemic models have been studied in the literature of mathematical epidemiology to describe the transmission dynamics of various infectious diseases such as measles, mumps, gonorrhea, West-Nile virus and HIV/AIDS ^[15,24,39].

The main purpose of this study is to estimate the transmission rate of HFMD among different age groups by using a multi-group model with population transfer and flow. The proposed model was analyzed by combining analytical and numerical techniques, focusing on the different types of HFMD case data in mainland China in 2014. The article is organized as follows. In the next section, we collect the surveillance data of HFMD in China and establish a multi-group HFMD model, then investigate the disease-free equilibrium and basic reproduction number of the system. The Section 3 presents the optimal parameters, simulation of the residentially-scattered children, childcare and student clinical infectious data from 2014. Sensitivity analysis of the basic reproduction number is carried out in Section 4. And we have conclusion and discussion in last section.

2.   The multi-group HFMD model

2.1.   Data

The Ministry of Health of the People's Republic of China declared that HFMD was ranked as a Category C Infectious Disease (Monitoring and Managing of Infectious Disease) on May 2nd, 2008. The Chinese Center for Disease Control and Prevention (China's CDC) collects confirmed case infected by HFMD which is in mainland China (i.e., except Hong Kong, Macao and Taiwan) ^[1,43] every month, there are teacher, farmer, nurse, medical staff, houseworker, residentially-scattered children, childcare, student and so on by occupation which divided into 19 classes. The case data from 2014 can be seen, total number of residentially-scattered children (73.58%), childcare (22.85%) and student (3.21%) was 99.64% ^[43]). There are data information includes the area code, gender, occupation, date of birth, address, date of onset, date of diagnosis, especially, classification of disease which is labeled as clinically diagnosed cases. The data were released and analyzed anonymously.

2.2.   Model formulation

The total population $N$ in our model is classified as $n$ groups: residentially-scattered children (the vast majority are less than 3 years old), childcare (about 3 to 6 years old), student (about 6 to 20 years old) and others (a large proportion of them are adults, almost over 20) as well as into five disjoint compartments for $n$ groups: susceptible $S$, exposed $E$ (infected but not infectious), clinical infectious $I$, subclinical infectious $L$ and recovered $R$. Every group moves from their susceptible compartments into the exposed compartments, where they display no symptoms and can not infect others. And there may be clinical or subclinical for the infectious individuals. After infection, all the individuals become recovered. Because of the incubation period (i.e., exposed state) and the duration of the HFMD (i.e., infectious state) were shorter, we only consider two categories compartments (i.e., susceptible and recovered state) has transfer of each other. That is to say, model (2.1) can express the residentially-scattered children go to school, the students enters a higher school, adults change jobs and the like.

Therefore, the following $n$-groups model is derived to describe the HFMD dynamics:

The parameters in the model are summarized in the following list:

● ${{\Lambda _i}}$: influx of individuals into the $i$-th group;

● ${\beta _{ij}}$: rate of disease transmission between susceptible individuals in group $i$ and clinical infectious individuals in group $j$;

● ${\gamma _{ij}}$: rate of disease transmission between susceptible individuals in group $i$ and subclinical infectious individuals in group $j$;

● ${\sigma _{ik}}$: transfer rate move from the $k$-th susceptible group into $i$-th susceptible group;

● ${\varepsilon _{il}}$: transfer rate move out the $i$-th susceptible group into $l$-th susceptible group;

● ${\kappa _{ip}}$: transfer rate move from the $p$-th recovered group into $i$-th recovered group;

● ${\omega _{im}}$: transfer rate move out the $i$-th recovered group into $m$-th recovered group;

● ${{\lambda _i}}$: remove rate from recovered to susceptible in group $i$;

● ${{\alpha _i}}$: rate of progression to infectious in group $i$;

● ${{\rho _i}}$: proportion of infective becoming clinical infectious in group $i$;

● ${{\delta _i}}$: recovery rate of clinical infectious individuals in group $i$;

● ${{\eta _i}}$: recovery rate of subclinical infectious individuals in group $i$;

● ${{\mu _i}}$: natural death rate in group $i$.

According to biological significance, parameters ${{\Lambda _i}}$, ${\sigma _{ik}}$, ${\varepsilon _{il}}$, ${\kappa _{ip}}$, ${\omega _{im}}, i, k, l, p, m = 1, 2, \cdots, n$, are all non-negative, and other parameters are all positive. The disease-free equilibrium (DFE) of model (2.1) is ${P^0} = \left({S_1^0, 0, 0, 0, 0, S_2^0, 0, 0, 0, 0, \cdots, S_n^0, 0, 0, 0, 0} \right) \in {R^{5n}}$, and ${S_i^0, i = 1, 2, \cdots, n}$ satisfied with the next algebraic equations.

Following P.V.D. Driessche and J. Watmough ^[27], we can compute the basic reproduction number. Note that the basic reproduction number $\mathcal{R}_0$ stands for the number of infected during the initial patient's infectious (not sick) period. $\mathcal{R}_0$ is used to determine whether a disease die out (if ${\mathcal{R}_0 < 1}$) or become epidemic (if ${\mathcal{R}_0 > 1}$), but for models with complex dynamics, ${\mathcal{R}_0 < 1}$ is not the only condition to guarantee that the disease is extinct, however the smaller the better ^[25,26]. The next-generation matrix approach in ^[27] was applied to calculate the basic reproduction number, $\mathcal{R}_0$. Rewriting the middle 3n equations of system (2.1) as $\dot x = \mathscr{F} - \mathscr{V}$, and $x = {({E_1}, {E_2}, \cdots, {E_n}, {I_1}, {I_2}, \cdots, {I_n}, {L_1}, {L_2}, \cdots, {L_n})^T}\in {R^{3n}}$. For this purpose, we can write the right-hand side of model (2.1) as $\mathscr{F - V}$ with

$\mathscr{F}$ and $\mathscr{V}$ are $3n$ dimensional column vectors. Denote $O$ is an $n$ dimensional null matrix. Calculating the Jacobian matrices, $\bf{F}$ and $\bf{V}$, at the DFE, we have

where,

where,

The basic reproduction number $\mathcal{R}_0$ is the spectral radius of $\bf{F{V^{ - 1}}}$.

3.   Parameter estimation and hypothesis

According to the data from China's CDC ^[1,43], More than 99% of clinical infectious $I$ are residentially-scattered children (young children who are not in school), childcare and student (elementary school students, middle school students, college students and so on). Choosing $n = 4$, then model (2.1) becomes a 4-groups model. Hence, the infectious class in our model is divided into four age groups: residentially-scattered children (${I_1}$), childcare (${I_2}$), student (${I_3}$) and others (${I_4}$, most of them are adults) and we'll just fit the data of ${I_1}, {I_2}, {I_3}$. Considering the actual situation, when residentially-scattered children arrive the school-age, they will go to kindergarten or nursery, then become childcare. Typically a few years later, they will go to primary school and become students. College students find jobs after they graduate. So we only estimate the transfer rates (essentially enrolment rates): ${\sigma}_{21}, {\sigma}_{32}, {\sigma}_{43}, {\varepsilon }_{12}, {\varepsilon }_{23}, {\varepsilon }_{34}, {\kappa }_{21}, {\kappa }_{32}, {\kappa }_{43}, {\omega }_{12}, {\omega }_{23}, {\omega }_{34}$, other transfer rates equal to 0.

By generating re-sample, we can product a larger artificial data that is generated based on the existing limited monthly data. Using the linspace function in Matlab, we interpolated the 12-month data and turned into 365-day data. In order to keep the total number of data, the interpolation formula as following:

where, ${{D_1}({s_i})}, i = 1, 2 \cdots, 12,$ denote the 12-month actual data, ${{D_2}({t_j})}, j = 1, 2 \cdots, 365,$ denote the 365-day data after the interpolation. ${\hat D_2}({t_j}), j = 1, 2 \cdots, 365,$ denote the 365-day data after conversion. With the aid of linear interpolation, we will obtain more useful data, and the fit results will be better.

Referring to literature ^[36], we select $\Lambda _1 = 4000$, $\Lambda _2 = \Lambda _3 = \Lambda _4 = 0$ (certainly, new babies are residentially-scattered children), ${\alpha _1} = {\alpha _2} = {\alpha _3} = 1/4.38, {\alpha _4} = 1/2.2$. One assume that the person's natural death follows a uniform distribution, then natural death rate is calculated as ${\mu _i} = 1/(74.83\times365), i = 1, 2, 3, 4,$ since life expectancy is 74.83 years old in China in 2014 ^[42]. The remaining 54 parameters in model (2.1) and 20 initial values are estimated which using the hybrid optimization algorithm by DEDiscover (an optimization software) ^[44], which has a superior performance over many alternative methods ^[29] and has been used in several previous studies ^{[30,31,32,33]}. For more details about DESQP, see Liang et al. ^[32]. The parameters and initial values of model (2.1) are listed in Table 2 and 3 in Appendix A. Moreover, the basic reproduction number was estimated as $\mathcal{R}_0 = 1.0328$ on the basis of our parameter values. Note that many previous literatures estimate the basic reproduction number in mainland China are between $1.0414$ to $1.7420$ ^{[34,35,28,37,38,36]} (seeing Table 1), and our estimate of the basic reproduction number is smaller than these results. By model (2.1), one carries on the data fitting to the number of clinical infectious (i.e., ${I_1}, {I_2}, {I_3},$), as shown in Figures 1, 2 and 3, the numerical results are found to be a good match with the data of HFMD in China in 2014.

We only fit the clinical data of three classes of patients according to model (2.1). The actual patient data has a double peak phenomenon, which we suspect is related to factors such as the opening time of school and vacation, the frequent flow of people, temperature and humidity, etc. If we want to characterize the double peak phenomenon, we usually need to use the dynamic model with periodic solutions, which will be an important issue for our future research.

4.   Sensitivity analyses

To examine the influence of parameter changes (the ranges given in Table 2 and 3), especially the rate of disease transmission (${\beta _{ij}}, {\gamma _{ij}}, i, j = 1, 2, 3, 4$), durations of latency before the onset of symptoms (${1 / {{\alpha _i}}}, i = 1, 2, 3, 4$), the duration of infectious following symptoms onset (${1 / {{\delta _i}}}, {1 / {{\eta _i}}}, i = 1, 2, 3, 4$), quarantine ratio (${\rho _i}, i = 1, 2, 3, 4$), on the control of the basic reproduction number $\mathcal{R}_0$. Following the method in ^[40] method, we carried out sensitivity analyses. A small perturbation ${\delta _\lambda }$ to a parameter $\lambda$ and the corresponding rate of change in $\mathcal{R}_0$ as ${\delta _{\mathcal{R}_0}}$

and normalized sensitivity index ${{\varphi _\lambda }}$ is defined as

the sensitivity indices of $\mathcal{R}_0$ are shown in Table 2 and 3 in Appendix B. The greater the absolute value of sensitivity indices of $\mathcal{R}_0$, the more sensitive to $\mathcal{R}_0$.

From sensitivity analyses of $\mathcal{R}_0$, we can obtain many useful conclusions. ${\beta _{11}}$, ${\beta _{21}}$, ${\beta _{31}}$ and ${\beta _{41}}$ are the most sensitive parameters among ${\beta _{1j}}$, ${\beta _{2j}}$, ${\beta _{3j}}$, ${\beta _{4j}}$, $j = 1, 2, 3, 4$ for $\mathcal{R}_0$, respectively. ${\gamma _{14}}$, ${\gamma _{24}}$, ${\gamma _{34}}$ and ${\gamma _{44}}$ are the most sensitive parameters among ${\gamma _{1j}}$, ${\gamma _{2j}}$, ${\gamma _{3j}}$, ${\gamma _{4j}}$, $j = 1, 2, 3, 4$ for $\mathcal{R}_0$, respectively. Apparently, ${\beta _{ij}}$ denotes the rate of disease transmission between susceptible individuals in group $i$ and clinical infectious individuals in group $j$, almost three quarters of clinical infectious are residentially-scattered children (${I_1}$); ${\gamma _{ij}}$ denotes the rate of disease transmission between susceptible individuals in group $i$ and subclinical infectious individuals in group $j$, and overwhelming majority of subclinical infectious individuals are adults (${I_4}$), this conclusion can also be referred to in the literature ^[28]. If we continue to compare ${\beta _{11}}$, ${\beta _{21}}$, ${\beta _{31}}$ and ${\beta _{41}}$, then ${\beta _{41}}$ is the most sensitive parameter. Similarly, ${\gamma _{44}}$ is the most sensitive parameter among ${\gamma _{14}}$, ${\gamma _{24}}$, ${\gamma _{34}}$ and ${\gamma _{44}}$. Whether the contact infection rate describes susceptible adults contact with clinical residentially-scattered sick children (${\beta _{41}}$), or describes susceptible adults contact with other subclinical adult patients ${\gamma _{44}}$, both illustrate the adults have a huge impact on the prevalence of HFMD. Not only that, from many previous literatures, the numerous people with subclinical adults infection who carry the HFMD virus but have no symptoms play an important part in leading to the pandemic ^[21,22,28]. On the basis of sensitivity indices of $\mathcal{R}_0$ (seeing Table 4 and 5), the parameters were divided into 4 categories. (1) The very sensitive parameters for $\mathcal{R}_0$ are: ${\gamma _{44}}, {\eta _4},$; (2) The more sensitive parameters for $\mathcal{R}_0$ are: ${\gamma _{14}}, {\rho _4}, {\delta _1}, {\beta _{41}}, {\rho _1}, {\delta _4}, {\gamma _{41}}, {\eta _1}$; (3) A little sensitive parameters for $\mathcal{R}_0$ are: ${\alpha _4}, {\beta _{21}}, {\gamma _{42}}, {\eta _2}, {\delta _2}, {\beta _{42}}, {\gamma _{43}}, {\eta _3}, {\rho _3}, {\gamma _{34}}, {\beta _{31}}$; (4) The remaining parameters are extremely insensitive.

5.   Conclusion and discussion

It is critical that every child is reached EV71 vaccine to ensure that population immunity is highly enough to prevent spread of HFMD virus and protect against further importations. But EV71 vaccine coverage is still very poor in China, we didn't consider vaccine for our model, so one won't discuss the EV71 vaccine in the following. On the basis of data analysis, we will put forward some control measures to reduce the spread of the epidemic outbreak.

(1) ${\gamma _{44}}$ is the most sensitive parameter for $\mathcal{R}_0$, in order to inhibit disease transmission one can reduce the contact rate between susceptible adults with subclinical adult patients. Hence, adults, especially parents, young teachers, nursing workers, medical staff and guardian of children should pay attention to personal hygiene. ${\eta _4}$ is the second sensitive parameter for $\mathcal{R}_0$, in order to effectively control the spread of HFMD can shorten the time of the adults germ-carrying (reduce $1/{\eta _4}$). Because sometimes adults are exposed to the virus but there are no obvious symptoms. As a matter of fact, there are 10 other parameters associated with adults, ${\gamma _{14}}, {\rho _4}, {\beta _{41}}, {\delta _4}, {\gamma _{41}}, {\alpha _4}, {\gamma _{42}}, {\beta _{42}}, {\gamma _{43}}, {\gamma _{34}}$, and they ranked the top 20 in the most sensitive parameters, which provides further evidence that adults are important in preventing the epidemic of HFMD.

(2) ${\delta _1}$ is a more sensitive parameter for $\mathcal{R}_0$, shorten treatment time for residentially-scattered sick children can inhibit the spread of HFMD. Children younger than three years old are high-risk group, accept agammaessive treatment (reduce $1/{\delta _1}$) for reducing highly pathogenic infection source and avoiding spreading again.

(3) ${\eta _1}$ is also a more sensitive parameter for $\mathcal{R}_0$, shortening the time of the residentially-scattered germ-carrying can effectively control the spread of HFMD. Therefore, it is an effective measure that we get great health-care education such as washing hands before meals and after using the toilet, and making air fresh indoors and so on.

(4) For childcare and student, we still have similar preventive and control measures. In contrast to adults, they are not key groups. Compared with residentially-scattered children, they are not high-risk groups.

This paper divides the patients with HFMD into four categories: residentially-scattered children (the vast majority are less than three years old), childcare, student and others (most of them are adults). The detailed data analysis shows that subclinical adults infection who carry the HFMD virus but have no symptoms play an important part in leading to the pandemic. It is consistent with the previous study ^[28]. This paper presents a more detailed classification of patients, and we confirmed that children under three years of age are indeed at high risk.

We conclude that the prevention of HFMD epidemic as follow. First, adults are the focus of attention, especially those who have more contacts with children. So they should be encouraged to pay a strict attention to personal hygiene, such as parents, young teachers, nursing staffs, medical staffs and child guardians in order to control the spread of HFMD effectively. Next, scattered children (under 3 years old) who are at high risk shall be given shorter treatment time to reduce the sources of highly pathogenic infections and to avoid retransmission. And then, for young children and students, they are not the key group, nor the most at-risk group, but the disease is still transmitted among them. Therefore, for everyone, we have a recommendation that we should strengthen health care education, deeply popularize hygiene and safety knowledge such as hand washing before meals, hand washing after going to the toilet, and keep the air fresh inside, so as to reduce the spread of HFMD.

Acknowledgments

We would like to thank anonymous reviewers for very helpful suggestions which improved greatly this manuscript. The work is supported by National Natural Science Foundation of China (Nos. 11471133 and 11547006), Chongqing Postdoctoral Science Foundation Special Funded Project (No. Xm2017139) and Scientific Research Project of Hubei Provincial Department of Education (No. B2017039).

Conflict of interest

The authors declare that they have no competing interests.

Appendix A

Appendix B