The impact of regular school closure on seasonal influenza epidemics: a data-driven spatial transmission model for Belgium

Background School closure is often considered as an option to mitigate influenza epidemics because of its potential to reduce transmission in children and then in the community. The policy is still however highly debated because of controversial evidence. Moreover, the specific mechanisms leading to mitigation are not clearly identified. Methods We introduced a stochastic spatial age-specific metapopulation model to assess the role of holiday-associated behavioral changes and how they affect seasonal influenza dynamics. The model is applied to Belgium, parameterized with country-specific data on social mixing and travel, and calibrated to the 2008/2009 influenza season. It includes behavioral changes occurring during weekend vs. weekday, and holiday vs. school-term. Several experimental scenarios are explored to identify the relevant social and behavioral mechanisms. Results Stochastic numerical simulations show that holidays considerably delay the peak of the season and mitigate its impact. Changes in mixing patterns are responsible for the observed effects, whereas changes in travel behavior do not alter the epidemic. Weekends are important in slowing down the season by periodically dampening transmission. Christmas holidays have the largest impact on the epidemic, however later school breaks may help in reducing the epidemic size, stressing the importance of considering the full calendar. An extension of the Christmas holiday of 1 week may further mitigate the epidemic. Conclusion Changes in the way individuals establish contacts during holidays are the key ingredient explaining the mitigating effect of regular school closure. Our findings highlight the need to quantify these changes in different demographic and epidemic contexts in order to provide accurate and reliable evaluations of closure effectiveness. They also suggest strategic policies in the distribution of holiday periods to minimize the epidemic impact. Electronic supplementary material The online version of this article (doi:10.1186/s12879-017-2934-3) contains supplementary material, which is available to authorized users.

We describe here in detail the inference procedure used to approximate fluxes of commuters per age class i = a, c in Belgium. For each commuting link l of the French commuting network, we computed the commuting distance d(l) and the fraction ρ(l) of commuters of age class i. We filtered all links having less than 30 commuters. We considered seven bins of distance, according to the definition used in the "Enquête National Transport et Déplacements 2008" (National survey on transport and mobility, 2008): [0, 2] km, ]2, 5] km, ]5, 10] km, ]10, 20] km, ]20, 40] km, ]40, 80] km, and > 80 km. The distribution obtained for each distance bin is then used to infer the fraction of Belgian commuters in age class i traveling on the same distance bin. A comparison between the empirical distributions obtained from the French commuting data and the reconstructed distributions for Belgium is shown in Figure S1. A good agreement is found for all distance bins, with a noisier behavior obtained for the bin class d(l) ≤ 2km, due to poor statistics. Plots also show that children commute at shorter distances than adults.  Figure S1: Probability distribution of the fraction of children commuters at specific distance bins: comparison between the empirical distributions obtained from the French commuting data and the reconstructed distributions for Belgium.

Contact Matrices
The average number of contacts made by participants in age class i = 1, 2 with people in age class j = 1, 2 is given by M ij . The per capita contact rates are then summarised in the contact rate matrix which is rescaled to a normalised contact matrix where N tot is the total population of Belgium We note here that the values C ij are scale invariant, that is We consider C ij to describe the social interaction in each patch, i.e. C (p) following Eq.(S1), with N (p) (t) being the total population in the patch at time t.

Details of the compartmental model in each patch
Each patch receives commuters from k in (p) patches and moves residents to k out (p) other patches, with k in and k out representing the indegree and outdegree, respectively, of patch p in the commuting network. Commuters are modeled with separate compartments, in order to track them in their movements from residence to destination and back. At each time step t, the population of a patch p is composed of the following subpopulations, each described by a two-age class SEIR disease progression model: • individuals who reside in patch p and do not commute: • k out (p) subpopulations of individuals who reside in patch p and commute to another patch q: , with q neighbor of p, and i = c, a; • k in (p) subpopulations of individuals who reside in a patch q and commute from patch q to patch p: , with q neighbor of p, and i = c, a.
Accounting for commuting, we can then write the force of infection for a susceptible individual of age class i in patch p and time t:

Influenza transmission
Here we describe influenza transmission for the subpopulation of residents of a given patch p (we drop the p for simplicity). The extension to the other subpopulations present in the patch is straightforward. The probabilities associated to SEIR transitions for age class i in a small enough time interval dt are given by: The number of individuals in age class i newly entering the E, I, and R class are extracted with binomial distributions (B):

Derivation of the next generation matrix
In calculating the values of R we disregard mobility. The in-patch model becomes therefore a two age-classes stochastic SEIR model. Its deterministic couterpart can be written as: Using Diekmann's approach we linearize the equations of the infectious compartments E, I around the disease free state with the correct immunity fraction a and obtain for the following system of linear equations restricted to the infectious compartments: The next generations matrix in the patch therefore reads (I I I is the identity matrix): which in components gives the result reported in the main text:

Computational details
The code of the simulations was written in C++ and made use of the Mersenne-trwister random generator and the binomial extraction procedures as provided by the Boost Libraries v1.58.0. Compiling was done using the gnu c++ compiler version 4.8.1 with optimization level 3. Table S1 presents the list of district names and associated IDs used in the study.

Calibration procedure
We minimized the Weighted Least Square function W LS(β n , α n ) computed on the median normalized incidence curves, considered from the start of the epidemic up to the peak time. The calibration is performed on Brussels district only and for each set of parameters (β n , α n ) we performed 1,000 simulations. Here β n is the explored per-contact transmission rate, and α n is the rescaling factor for the simulated incidence in Brussels district to account for possible sampling biases in the initial condition. The calibration is performed on normalized incidence curves to discount the effects of unknown GP consultation rates.
To reduce the number of points to explore and cope with stochastic fluctuations we considered iterative resampling through a particle filter/bootstrap method. For each level l we calculate a weight distribution for each (β n , α n ) as follows: which allows us to define a filtering/resampling transition probability: p(β, α; l + 1 | β l 1 , . . . , β l n , α l 1 , . . . , α l n ; l) = β l n ,α l n w l (β l n , α l n )V βn,αn) (β, α) (S13) where V (a0,b0) (a, b) is the uniform distribution over the Voronoi cell centred in (a 0 , b 0 ). We can then resample N -particles at level l + 1 given the M -particles at level l and repeat the process iteratively until the filtering Ath Aat Table S1: List of district names and associated IDs.
probability is almost uniform, and therefore the filter does not work any more. Here we used 20 particles at each level (except the first) and then we stopped when the number of effective particles defined as: (S14) was greater of 19, which corresponds to uniformity of filtering probability. For each level l, we thus obtain a set of 20 pairs (β l n , α l n ) whose distribution is used to estimate the values of α and β that minimize W LS(β, α).

Additional validation results
Calibration results are listed in Table S2. Figure S2 shows the comparison in the peak timing between simulations calibrated with values of Table S2 and surveillance data.   Figure S2: Left: Boxplot of the peak time difference ∆T d per district between simulations and empirical data. Numbers represent Belgium districts, see Table S1 for corresponding names. Right: Geographical map of the median peak time difference per district.