Threshold control strategy for a Filippov model with group defense of pests and a constant-rate release of natural enemies

1.   Introduction

Hand-foot-and-mouth disease (HFMD) causes a substantial disease burden in Asian regions including China, mainly in children below 5 years old and partly in adolescent students ^[1,2]. Enterovirus 71 (EV71), coxsackievirus A16 (CA16) are the most common enterovirus serotypes causing HFMD in China ^[3,4]. Although most cases of HFMD are mildly self-limiting, characterized by fever, flat spots or bumps on the hands, feet and mouth, possible transmission routes of HFMD contain casual contact with contaminated discharges, fluid of blisters or stool from infected persons, or contaminated objects, severe cases may develop into heart and lung failure and even death ^[5,6,7].

In the last decade, many severe HFMD outbreaks have occurred in East Asia and Southeast Asia, which has become a serious public health problem in the affected countries and regions. Despite generation and approval of an EV71 vaccine in China in 2016, recently several HFMD outbreaks were caused by other pathogens, making it important to take additional actions to control HFMD epidemics ^[8]. The increasing burden of HFMD has become a serious public health problem in China, causing the government to pay greater attention to the risk factors of this disease, and focus on the effective prevention and control strategies.

A distinct seasonal pattern with annual peaks in the summer was reported in temperate regions ^[9], whereas in tropical and subtropical regions, more than two peaks may occur within a year ^{[10,11,12,13,14]}. Many studies have confirmed that meteorological factors have an important influence on the incidence of HFMD despite the inconsistent results from different studies ^{[9,12,15,16,17,18,19]}. Relative-humidity and temperature are regarded as the top two meteorological factors affecting the HFMD transmission ^{[15,17,18,19]}. Attending kindergarten, child-care center or school and public gathering were also pointed out as risk factors for HFMD ^[20,21,22]. However, it is not clear whether above factors are related to multiple-peak pattern, and when and how these factors play a critical role in determining the HFMD transmission in one year.

Mathematical epidemiology is the study of the spread of disease in space and time, with the aim to identify the factors that are responsible for or contributing to their occurrence. The basic reproduction number ($R_{0}$) is often used to assess the transmission potential or transmissibility of an infectious disease, which indicates the average number of expected secondary infections resulting from a single infectious case ^[23]. In general, if it is greater than one, the disease will continue to spread, while if it is less than one, the disease will eventually disappear ^[24,25]. Both statistically robust analytical methods (model-free estimation) and dynamic modelling of infectious diseases (model-based estimation) have contributed greatly to estimate $R_{0}$, which may be used to guide infectious disease prevention, control and eradication ^[10,26,27].

Now, mathematical models have become one of the important tools for analyzing the spread and control of HFMD ^{[3,28,29,30,31,32,33]}. The aim of this study is to investigate the transmission dynamics and examine when and how its potential factors play a decisive role in the annual multiple-peak pattern based on the surveillance data of HFMD in Wenzhou city, Zhejiang province, China. We first estimated daily and annual model-free basic reproduction number $R_{0}^{free}(day)$ and $R_{0}^{free}(year)$ based on the infected data and the infection characteristics of HFMD, in which $R_{0}^{free}(day)$ means the average number of secondary cases that each infected individual would infect if the conditions remained as they were at the time of a day, and $R_{0}^{free}(year)$ is the annual average of $R_{0}^{free}(day)$. Then we developed a mathematical model incorporating the potential factors, which are screened and identified by the correlation with the obtained daily model-free basic reproduction number, including school opening or not, and meteorological factors. After then, reported HFMD infection data were further used to estimate the undetermined model parameters. Thus, the annual model-based basic reproduction number $R_{0}^{model}(year)$ can be calculated after some reasonable approximation, and sensitivity analysis of the model was performed to estimate the resultant effects of the potential factors. In the end, some feasible control strategies of HFMD were proposed based on the above analyses.

2.   Materials and methods

2.1.   Data collection

Infected data of HFMD in Wenzhou between 1 January 2010 to 31 December 2015 were collected by Wenzhou Center for Disease Control and Prevention ^[34], including basic demographic characteristics of HFMD cases, and daily incident cases. According to the obtained records, all the cases were children aged between 0 and 14 years old, and there were 203,488 cases in all.

School summer and winter vacation as well as school opening duration was determined according to the guidelines of the school calendar issued by the Education Bureau ^[35].

Daily meteorological factors, including minimum temperature, mean temperature, maximum temperature, minimum relative humidity and mean relative humidity, were collected from the China Meteorological Data Sharing Service System ^[36]. A weather monitoring station is located in Ruian County, Wenzhou City (N27^o', E120^o39^'), from which the data are usually used to represent the meteorological conditions of the whole Wenzhou city.

2.2.   Estimation of the daily model-free basic reproduction number

To quantify the intensity of transmission over time, we first suppose that serial interval (the time between the onset of symptoms in a primary case and the onset of symptoms in secondary cases) is equal to the sum of $L$ and $D$, where $L$ represents the duration of incubation and $D$ denotes the duration of infection ^[37]. Then serial interval distribution was obtained by calculating the joint probability density function of $L$ and $D$. Thus, according to the method in ^[26], we can estimate the daily model-free basic reproduction number $R_{0}^{free}(day)$ by using "epiEstim" Excel software package ^[26]. According to the annual average, we can obtain the annual model-free basic reproduction number $R_{0}^{free}(year)$. See Supplementary A for more details.

2.3.   Dynamic model of HFMD

According to the primary epidemiological properties of HFMD and the implemented public health interventions (quarantine, isolation and hygiene precaution), we divide the total population ($N$) into five compartments, named as susceptible ($S$), exposed but not yet infectious ($E$), infectious ($I$), quarantine ($Q$) and recovered ($R$). Figure 1 illustrates the flowchart of the five compartments, which describes how individuals can move among the states.

Based on the flowchart in Figure 1, the dynamic model can be described by the following non-autonomous differential equations:

The meaning and values of parameters, and initial conditions in model (2.1) are given in Table 4 (See Supplementary B). Note that once individuals are infected with the virus, even though they are growing in age, they can only leave the system after evolving into recovered compartment. Thus, the parameter of transfer rate due to age growth $c$ only appears in the first, fifth and last equations in (2.1).

Here $\beta(t)$ is the contact rate between susceptible ($S$) and infectious ($I$), which is a time-varying, continuous, positive function. Concretely, $\beta(t)$ can be divided into five sections: (1) the interaction factor between school opening and meteorological effect during spring semester, $\beta_{s\times m}(t)$; (2) the interaction factor between school opening and meteorological effect during fall semester, $\beta_{f\times m}(t)$; (3) the meteorological factor during spring semester, $\beta_{ms}(t)$; (4) the meteorological factor during fall semester, $\beta_{mf}(t)$; (5) the meteorological factor during winter and summer vacation, $\beta_{mr}(t)$. Therefore, $\beta(t)$ can be expressed as

For detailed expressions of (2.2) and more details, see Supplementary B.

2.4.   Simulations

Based on the reported HFMD cases and the dynamic model (2.1), an adaptive Metropolis-Hastings algorithm was conducted to carry out the Markov Chain Monte Carlo (MCMC) procedure to estimate the undetermined parameters in contact rate and their standard deviations (see Supplementary C for details). The algorithm runs for $1000000$ iterations with a burn-in of $500000$ iterations, with the Geweke convergence diagnostic method employed to assess the convergence ^[38]. The mean values of the estimated parameters and their standard deviations are listed in Table 5. Note that the fluctuations of meteorological factors in different years can be approximately regarded as periodic changes. i.e., the contact rate $\beta(t)$ can be approximately regarded as annual periodic function. Thus, based on the results in ^[39,40], the annual model-based basic reproduction number $R_{0}^{model}(year)$ can be calculated by estimated parameters and simulation numerical integration (See Supplementary D for details).

2.5.   Sensitivity analysis

Note that there are $19$ undetermined parameters in contact rate, each with different mean value and variance, which means that it is tedious and extremely complex to study their influence on the outcomes. In order to detect the underlying factors more simply, a better idea is to do the sensitivity analysis for different categories of parameters. Thus, according to the expression (2.2) of contact rate, we classified all parameters into five factors, also marked as $\beta_{s\times m}$, $\beta_{f\times m}$, $\beta_{ms}$, $\beta_{mf}$ and $\beta_{mr}$ (see Supplementary B for details). Concretely, we use Latin hypercube sampling (LHS) and partial rank correlation coefficients (PRCCs) ^[41,42] to study the dependence of the above five factors on daily contact rate and annual model-based basic reproduction number and daily contact rate from 2010 to 2015.

Firstly, to identify key underlying factors that influenced the HFMD dynamics, we used LHS and PRCCs to examine the sensitivity of the annual model-based reproduction number to the five factors ^[42,43]. Normal distribution for all input parameters with the mean values and standard deviations was given in Table 5. By substituting the sampled parameters' values into the formulas (7.2-7.6) in Supplementary B, all five factors' values were obtained. Then the annual global PRCCs were calculated between the five factors and the corresponding output variable, i.e., the annual model-based reproduction number.

Secondly, because the correlation between daily model-free basic reproduction number and daily contact rate has statistical significance, time-varying sensitivity was done between daily contact rate and five factors. Concretely, we still choose all five factors as the input variables, but use the estimated daily contact rate as the output variable. Using similar sampling parameter method as above, time-varying PRCCs were calculated for each factor. Here the sample size $N_{s}$ is $10000$ in both global and time-varying sensitivity analyses. Similar to ^[41,42], we considered absolute values of PRCC $> 0.4$ as indicating an important or strong correlation between input parameters and output variables, values between $0.2$ and $0.4$ as moderate or weak correlations, and values between $0$ and $0.2$ as not significantly different from zero.

2.6.   Statistical analysis

For non-normal distribution data, such as all meteorological factors, median with the corresponding first and third quartiles (Q1-Q3) was reported. Otherwise, Mean $\pm$ standard deviation (SD) was presented, such as annual model-free and model-based basic reproduction number. Spearman correlation coefficient was used to measure the correlation between meteorological factors, and the correlation between daily model-free basic reproduction number and contact rate. According to whether the daily model-free basic regeneration number is less than one, we divide the meteorological factors into two groups ($R_{0}^{free}(day)\geq 1$ and $R_{0}^{free}(day) < 1$), and Mann-Whitney test was used to analyze the difference between the two groups. Furthermore, independent-samples t test was performed to analyze the difference between annual model-free and model-based basic reproduction numbers. In order to compare with the threshold value (one), one-sample t test was performed for both annual model-free and model-based basic reproduction numbers. A two-tail $P < 0.05$ was considered statistically significant.

3.   Results

3.1.   Daily model-free basic reproduction number and multiple-peak pattern of HFMD outbreaks

Figure 2 showed the estimated daily model-free basic reproduction number ($R_{0}^{free}(day)$) and the reported daily HFMD infectious cases in Wenzhou between 1 January 2010 to 31 December 2015. In every year, we can see that the basic reproduction number of HFMD (blue solid line in Figure 2) began to increase above unit during early March and it usually remained above one for around four months from early March to late June. After that, a rebound of basic reproduction number occurred between September and October in autumn, then after wards fluctuated below one for most of the time during winter. This is consistent with the multiple-peak pattern of HFMD outbreaks, which can be seen from the reported data (orange and green histogram in Figure 2).

3.2.   Daily model-free basic reproduction number and school opening

Figure 3 revealed that daily model-free basic reproduction number changed with school opening days during spring semester and fall semester, respectively. In the early stage of the spring semester, the basic reproduction number began to rise exponentially and quickly exceeded threshold one in no more than half a month (Figure 3A), and it usually peaked around the 30th day after school beginning and then remained above one until about the 120th day. Similarly, at the beginning of the fall semester, the basic reproduction number also grew rapidly and exceeded one in half a month (Figure 3B), and it usually peaked around the 30th day after the beginning of the fall semester and then fell down to less than one in a month. After two months of school beginning, the basic reproduction number began to increase again and reached a peak slightly above one, and then continued to decline through winter.

3.3.   Meteorological factors and daily model-free basic reproduction number

Various meteorological factors were found to be inter-related with each other. Table 1 summarized the correlation coefficients between all these factors. It was found that mean temperature, maximum temperature and minimum temperature were highly correlated with each other, with the correlation coefficients larger than $0.90$. Moreover, mean relative humidity and minimum relative humidity were also highly interrelated with a correlation coefficient $0.910$.

Table 2 showed the comparison of each meteorological factors between the group of daily model-free reproduction number $R_{0}^{free}(day)\geq 1$ and $R_{0}^{free}(day) < 1$. There were significant differences in all factors between the two groups ($P < 0.05$). Hence, combining the previous correlation analysis results, we included only mean temperature and mean relative humidity in the dynamic model, and excluded other highly correlated variables.

3.4.   Fitting the reported HFMD data and estimating the annual model-based basic reproduction number

We employed an adaptive Metropolis-Hastings algorithm to carry out extensive MCMC simulations ^[38], and to estimate mean values of parameters in the contact rate expression (2.2). Using large sample realization, we fitted the model (2.1) to the data of reported HFMD cases and the fitted parameters are shown in Table 5 in Supplementary C.

Figure 4 illustrates the best-fit numerical simulation of the model with the fitted parameters, which shows that the model was well fitted to the reported cases, and can provide the overall trend of HFMD infected population in Wenzhou city. Furthermore, the annual model-based basic reproduction number $R_{0}^{model}(year)$ from 2010 to 2015 were given in Table 3. Comparing with the threshold one, we find that all annual model-based and model-free basic reproduction numbers were significantly larger than threshold one, which indicates the consistency of two estimation methods. This suggests that our findings may be a possible mirror of the real transmission of HFMD i.e., HFMD is an endemic disease in Wenzhou city, which is consistent with the relevant studies in other regions of China ^{[29,30,31,33]}. However, we can also find that the model-based estimates are usually less than model-free estimates (Table 3), which should arouse our attention and further research in the future.

3.5.   Sensitivity analysis of the annual model-based basic reproduction number

Taking annual model-based reproduction number as the target output variable, we calculated the PRCCs of all five factors ($\beta_{s\times m}$, $\beta_{f\times m}$, $\beta_{ms}$, $\beta_{mf}$ and $\beta_{mr}$) in contact rate. From Figure 5, we found that all five factors were sensitive to the annual model-based basic reproduction number every year. As expected, the stronger the factor, the greater annual basic reproduction number. In particular, we find that the PRCCs of each factor change little every year, which indicates that each factor has better robustness to the impact of HFMD outbreak. By comparing PRCCs of each factor, we can further conclude that the influence of factors ($\beta_{s\times m}$, $\beta_{f\times m}$ and $\beta_{ms}$) are greater than that of factors ($\beta_{mf}$ and $\beta_{mr}$). These may provide potential guidance for HFMD prevention and control in strategy formulation.

3.6.   Time-varying sensitivity analysis

Firstly, the correlation between the daily model-free basic reproduction number and the daily model-based contact rate was demonstrated. Figure 6A shows the time series diagram of estimated daily model-free basic reproduction number and model-based contact rate. Furthermore, the results of Spearman correlation analysis indicates that there is a significant correlation between them ($\rho = 0.446$, $P < 0.05$).

Due to the positive correlation between them, time-varying analysis is performed through the time-varying PRCCs between the above-mentioned five factors and daily contact rate $\beta(t)$ with time (Figure 6B). Because the PRCCs of meteorological factor ($\beta_{mr}$) is very large during summer vacation, we can naturally conclude that the cases during this period are mainly caused by the effect of meteorological factors. In addition, during the spring semester, the PRCCs of meteorological factor ($\beta_{ms}$) is increasing gradually, but the PRCCs of interaction factor between school opening and meteorological effect ($\beta_{s\times m}$) is decreasing gradually, and those PRCCs almost maintain a high level throughout the spring semester. This may be the reason for the observed persistent peak of HFMD outbreaks in spring semester. In autumn semester, the variation trend of PRCCs of $\beta_{mf}$ and $\beta_{f\times m}$ is similar to $\beta_{ms}$ and $\beta_{s\times m}$, but the time for the interaction factor to maintain a higher level is obviously shorter than that for the meteorological factors to maintain a higher level. This may be the reason for the observed multiple small peaks in autumn semester. As a result, these may provide a quantitative characterization of the underlying factors in the annual multi-peak pattern of HFMD outbreak.

4.   Conclusions

Similar to the pattern in other subtropical and tropical regions, HFMD epidemics in Wenzhou occurred at least twice a year. The underlying factors for the annual pattern of HFMD multiple outbreaks have been investigated widely ^{[10,11,12,13,14]}. Some studies suggested a strong association between HFMD and weather ^[11,12]. In Hong Kong, researchers found that HFMD transmission was mainly associated with depletion of susceptible and limited impact were suggested from meteorological factors and school holidays ^[44]. However, in our study, it is clear that high contact rate was synchronized with the rapid increase of the daily model-free basic reproduction number, especially at the beginning of the school semester (Figure 6A). This pattern may be driven by a new cohort of susceptible children and few infected students entering kindergartens as well as the increase in contact rates among children in the new school semester ^[44]. Nevertheless, further confirmatory investigation is needed to demonstrate the impact of new students in kindergartens on HFMD transmission.

Though many studies indicated that meteorological factors play an important role in HFMD transmission ^{[9,12,15,16,17,18,19]}, this study further find that there are different impact effect in different seasons and their interaction with school opening. Meteorological factors in spring semester, the interaction factor between school opening and meteorological effect, made more significant contributions to HFMD transmission (Figure 5 and 6). Meanwhile, both daily model-free basic reproduction number and contact rate in summer holidays were obviously higher than that in winter holidays (Figure 6A). Those results were consistent with the previous studies that temperature and humidity may jointly affect this childhood disease ^[45], and such interactive impact needs to be considered when evaluating at different school semesters.

Based on the whole analysis, we can identify two feasible strategies for efficiently lessening the basic reproduction number, i.e., for containing HFMD epidemics. The first measure is to increase social distance. The transmission of enterovirus and coxsackievirus is most efficient in crowded settings and, therefore, most countries and regions, including Malaysia, Singapore, Taiwan, Hong Kong, and China, have adopted social distancing measures, such as closures of childcare facilities and schools, and cancellation of public functions involving children ^[46,47]. In our study, the optimum timing for implementation was at the beginning of every school semester when the peer-group at school started to gather together. At this time, as soon as an infected was detected in the childcare facilities and schools, social distancing measures must be immediately implemented. Such control measure may lower the basic reproduction number and decrease the peak HFMD incidence. Secondly, health education focusing on personal hygiene and good sanitation, including frequent hand washing, proper disposal of soiled nappies, and disinfection of soiled surfaces with chlorinated (bleach) disinfectants, should be addressed in the spring semester. There are two possible reasons why we highlight this measure. On the one hand, physical activity of adolescence is lower during winter and increases during warmer months and therefore, weather conditions may be associated with behavioral patterns that lead to increased contact among children, thereby facilitating the spread of HFMD infection ^[9,48,49]. On the other hand, weather changes have been credited with improving the duration of viral survival and virulence in the environment ^[50,51,52], and such changes increase the number of opportunities for hosts to become infected as a result of contact with contaminated surfaces, aspiration of contaminated fomites, or inhalation of aerosols produced by the coughs of infected individuals ^[53,54]. More evidence on the underlying biological mechanism of HFMD transmission is needed to improve the understanding of the potential effect of meteorological factors.

A few limitations should be considered when interpreting the findings from this study. First, the estimation of model-free basic reproduction number was based on hospitalized HFMD cases and may not reflect the trend of population incidence reliably. This may be the reason why model-free annual basic reproduction number had a minor difference with model-based annual reproduction number. Second, we adopted the dates of summer and winter vacation from the Education Bureau, but some kindergartens started schools earlier. This may partly explain the rebound of model-free basic reproduction number slightly earlier than the start of the school year in September. Finally, despite that we developed a mathematical model to establish causation between HFMD transmission, we may not incorporate all potential driving factors.

Acknowledgments

The authors are very grateful to the anonymous referees for their valuable comments and suggestions. This research was supported by the National Natural Science Foundation of P. R. China (Grant Nos. 11771448, 61672013, 61772017) and Huaian Key Laboratory for Infectious Diseases Control and Prevention (HAP201704).

Conflict of interest

The authors declare that they have no competing interests.

Supplementary

Supplementary A: Estimation of the daily and annual model-free basic reproduction number

Similarly to ^[37], we assumed that serial interval (the time between the onset of symptoms in a primary case and the onset of symptoms in secondary cases) is equal to the sum of $L$ and $D$, where $L$ represents the duration of incubation and $D$ denotes the duration of infection. As a result, ^[55] obtained the probability density function of incubation period of HFMD whose best-fitted function was gamma distribution. Concretely, since the estimated median incubation period was 5.4(95% CI 4.4-6.5) days, based on the parameters in ^[55], we know the probability density function of incubation $L$ was $f_{Incubation}(x) = \frac{b_{1}^{a_{1}}}{\Gamma(a_{1})}x^{a_{1}-1}e^{-b_{1}x}$ if $x\geq 0$, and $f_{Incubation}(x) = 0$ if $x < 0$, in which $a_{1} = 6.9637, b_{1} = 1.2217$.

In addition, the disability length for symptomatic HFMD or infection was estimated at 5.7(95% CI 5.4-6.0) in ^[56,57]. Suppose that HFMD infection is also obeys a gamma distribution. Then, based on the parameters in ^[56,57], the probability density function of infection $D$ was $f_{Infection}(x) = \frac{b_{2}^{a_{2}}}{\Gamma(a_{2})}x^{a_{2}-1}e^{-b_{2}x}$ if $x\geq 0$, and $f_{Infection}(x) = 0$ if $x < 0$, in which $a_{2} = 1.0264, b_{2} = 0.1901$.

Hence the joint probability density function of the sum of incubation period and infection duration ($L+D$), described as the probability density function of serial interval ^[26], was given by

Due to the difficulty in solving the above integral, we used numerical integration method to obtain the result (see Figure 7). After that, the daily model-free basic reproduction number ($R_{0}^{free}(day)$) can be calculated based on the infected data and the probability density distribution function of serial interval through the method and Ecxel software package "epiEstim" in ^[26]. Taking the annual average of $R_{0}^{free}(day)$, we can obtain the annual model-free basic reproduction number $R_{0}^{free}(year)$.

Supplementary B: Model parameters, initial conditions and contact rate

B1. Model parameters and initial conditions in model (2.1)

According to Statistical Year book and Census data from Wenzhou Statistical Bureau ^[58], the average daily number of new babies from 2010 to 2015 is $\lambda = 334.459$, and average daily mortality rate is $d = 1.904\times 10^{-5}$. Assume that the average daily deaths of people under 14 years old is $d\cdot N(t)$, the average daily number of people transferred out due to age growth is $c\cdot N(t)$, we can get $\frac{d N}{d t} = \lambda-(d+c)N$. By solving above equation, we obtained that $N(t) = \frac{\lambda}{d+c}+(N(0)-\frac{\lambda}{d+c})e^{-(d+c)t}$, and $N(0)$ is the initial susceptible. Fortunately, according to the Statistical Year book and Census data from Wenzhou Statistical Bureau, we find that the total population under 14 years in January 1, 2011 and December 31, 2015 is $1305373$ and $1355698$, respectively. Using these information, we can estimate the values of the parameters with the help of least square method. As a result, $c = 2.113\times 10^{-4}$ and $N(0) = 1292553$. Because the number of susceptible is almost equal to total population, we have $S(0) = N(0)$.

Overall, the definitions and values of all parameters and initial values of model (2.1) are listed in Table 4.

B2. Contact rate

The contact rate between susceptible and infected, $\beta(t)$, is a time-varying, continuous, positive function, which is mainly determined by school opening or not, and meteorological factors. The former can be divided into three phases, spring semester, fall semester and vacation (summer and winter vacation), while the mean temperature and relative humidity were selected to represent meteorological factors. For the sake of narrative convenience, several variables are expressed as bellow: $\delta_{r}(t)$ indicates whether it is the school vacation or not ($1$ represents school vacation while $0$ denotes school opening); $\delta_{s}(t)$ indicates whether it is the spring semester or not ($1$ represents spring semester (including the weekend and other short-term holiday during spring semester) while $0$ means not); $\delta_{f}(t)$ indicates whether it is the fall semester or not ($1$ represents fall semester (including the weekend and other short-term holiday during fall semester) while $0$ means not); $x_{s}(t)$ refers to the cumulative opening days of spring semester in one year; $x_{f}(t)$ refers to the cumulative opening days of fall semester in one year. Since we can see from Figure 3 that the both variation between daily model-free basic reproduction number $R_{0}^{free}(day)$ and school opening days in spring and fall semesters are similar to the form of probability density function similar of gamma distribution, the effect of the cumulative school opening days on the daily contact rate is defined as $y_{s}(t) = x_{s}(t)^{b_{1}-1}e^{-b_{2}x_{s}(t)}$ in spring semester and $y_{f}(t) = x_{f}(t)^{b_{3}-1}e^{-b_{4}x_{f}(t)}$ in fall semester.

Since ^[59,60] points out that there is no time lag effect between temperature and disease outbreak, to normalize the mean temperature per day, $T(t)$ is defined as the average temperature per day divided by $40$ because the highest daily mean temperature of Wenzhou city is no more than $40^{o}C$. However, ^[59] indicates that there exists time lag effect between humidity and disease outbreaks. Note that the daily model-free basic reproduction number $R_{0}^{free}(day)$ is a measure of the severity of HFMD outbreaks. We first calculate the correlation coefficient $\rho_{\bar{H}_{R}(t-\tau), R_{0}^{free}(day)}$ between $R_{0}^{free}(day)$ and the measured mean relative humidity with different time lag days $\bar{H}_{R}(t-\tau)$, $\tau = 0, 1, \cdot\cdot\cdot, 30$ (Figure 8). Then, according to the magnitude and stationarity of the correlation coefficients, we get that the outbreak of disease at time $t$ is mainly driven by the average relative humidity between $t-12$ and $t-23$ days. Thus, the average relative humidity at time is defined as $H_{R}(t) = \sum\limits_{\tau = 12}\limits^{23}r(\tau)\bar{H}_{R}(t-\tau)$, in which $r(\tau) = \rho_{\bar{H}_{R}(t-\tau), R_{0}^{free}(day)}/\sum\limits_{s = 12}\limits^{23}\rho_{\bar{H}_{R}(t-s), R_{0}^{free}(day)}$ is corresponding weight coefficient.

In summary, according to the seasonal variation, school opening and their interaction effect, we divide contact rate $\beta(t)$ approximately into the following five parts:

(1) the interaction factor between school opening and meteorological effect during spring semester,

(2) the interaction factor between school opening and meteorological effect during fall semester,

(3) the meteorological factor during spring semester,

(4) the meteorological factor during fall semester,

(5) the meteorological factor during winter and summer holiday,

Therefore, $\beta(t)$ can be expressed as (2.2) in the main text.

Supplementary C: Parameter estimation in contact rate

The variance of measured components, $I(t)$, was given by an inverse gamma distribution with hyper-parameters $(0.01, 4)$, where $0.01$ is the initial error variance which is updated by the inverse gamma distribution ^[61], and the small MCMC package provided in this website was used to estimate the parameters. When estimating unknown parameters and initial conditions for model (2.1), the following prior information were given: $b_{i}\in [0, +\infty), i = 1, 2, 3, 4$ and $b_{j}\in ( = \infty, +\infty), j = 5, 6, \cdot\cdot\cdot, 19$; the proposed density was chosen to be a multivariate normal distribution. These ranges were used to ensure good convergence of the MCMC chain.

By fitting model (2.1) to the reported HFMD cases from 2010 to 2015 in Wenzhou city, we estimate mean values of $b_{i}, i = 1, 2, \cdot\cdot\cdot, 19$, and their standard deviations (SD) in Table 5. Figure 9 is the distribution diagram of each parameter.

Supplementary D: Annual model-based basic reproduction number

According to ^[40], the annual model-based basic reproduction number $R_{0}^{model}(year)$ of model (2.1) can be defined by the spectral radius of the next infection operator. Note that the disease-free steady state of model (2.1) is $E_{0} = (\lambda/(d+c), 0, 0, 0, 0, \lambda/(d+c))$, the contact rate $\beta(t)$ can be approximated as an annual periodic function and the period is marked as $\omega$. Using the notation and definition in ^[40], as well as the general calculation procedure in ^[39], for the approximate periodic model (2.1) we have

In order to characterize model-based basic reproduction number of model (2.1), we consider the following linear ω-periodic system

with parameter $\theta\in (0, \infty)$. Let $W(t, s, \theta), \; t\geq s, \; s\in \mathbb{R}$, be the evolution operator of the system (7.7) on $\mathbb{R}^{3}$. For each $\theta\in (0, \infty)$, the matrix $-V(t)+\frac{F(t)}{\theta}$ is cooperative. Thus, the linear operator $W(t, s, \theta)$ is positive in $\mathbb{R}^{3}$ for each $t\geq s, \; s\in \mathbb{R}$. By Theorem 2.1 in ^[40], $R_{0}^{model}(year)$ can be defined as the unique solution of $\rho(W(\omega, 0, \theta)) = 1$ and $R_{0}^{model}(year) = \theta$ in this case, in which $\rho(W(\omega, 0, \theta))$ is the spectral radius of $W(\omega, 0, \theta)$.