Model-Based Comprehensive Analysis of School Closure Policies for Mitigating Influenza Epidemics and Pandemics

School closure policies are among the non-pharmaceutical measures taken into consideration to mitigate influenza epidemics and pandemics spread. However, a systematic review of the effectiveness of alternative closure policies has yet to emerge. Here we perform a model-based analysis of four types of school closure, ranging from the nationwide closure of all schools at the same time to reactive gradual closure, starting from class-by-class, then grades and finally the whole school. We consider policies based on triggers that are feasible to monitor, such as school absenteeism and national ILI surveillance system. We found that, under specific constraints on the average number of weeks lost per student, reactive school-by-school, gradual, and county-wide closure give comparable outcomes in terms of optimal infection attack rate reduction, peak incidence reduction or peak delay. Optimal implementations generally require short closures of one week each; this duration is long enough to break the transmission chain without leading to unnecessarily long periods of class interruption. Moreover, we found that gradual and county closures may be slightly more easily applicable in practice as they are less sensitive to the value of the excess absenteeism threshold triggering the start of the intervention. These findings suggest that policy makers could consider school closure policies more diffusely as response strategy to influenza epidemics and pandemics, and the fact that some countries already have some experience of gradual or regional closures for seasonal influenza outbreaks demonstrates that logistic and feasibility challenges of school closure strategies can be to some extent overcome.


Model description
We performed our analysis focusing on the United Kingdom by using an individual-based stochastic simulation model similar to the one developed for Europe in [1][2][3]: highly detailed country-specific socio-demographic data on age structure, household size and composition, employment rates by age, educational system were used to generate a synthetic population where each individual is given an age and is associated to a specific household and an occupation chosen between student, worker and inactive (unemployed/retired). The epidemic model is a spatially explicit, discrete time SEIR model with explicit transmission in households, schools and workplaces and in the general community. Throughout the whole course of the simulation cases from abroad were imported in the UK on the basis of the number of daily airport arrivals from the rest of the world, similarly to [1,4,5].
In the present version, the implemented school organization reflects the real situation existing in the UK in terms of both partitioning of students into different educational levels and class and school size. More in detail, four educational stages have been considered, following the educational structure in England and Wales, and quite in agreement with the situation of Scotland and Northern Ireland: pre-primary school (kindergartens, day-care centers) up to 5 years old, primary school for children from 5 to 11 years old, secondary school from 11 to 18 years old, higher education (post-secondary training, university, doctoral programs) for students older than 18. For each school level the corresponding information on school size has been taken into account in order to obtain a realistic size distribution [6]; in addition, primary and secondary schools have been further partitioned into classes and grades, reflecting information on average size by school level in the UK, namely 25 pupils per class in primary education and 20 students per class in secondary education [7]. All this information have been included in the model in order to build schools of the correct size and to assign every student to the correct school level (from pre-primary to higher education) and grade according to age; moreover, in order to apply mitigation policies based on gradual closure of schools, for primary and secondary schools students belonging to the same school were sub-grouped into classmates, same grade but different class or different grade.

Transmission model
As in [2][3][4]8] interaction between individuals (and thus infection transmission) can occur in four different settings, namely household, school, workplace or general community. Here we use different, setting-specific infectivity profiles over time and generation times. Moreover, as a consequence of the refinement of primary and secondary schools accounting for classes and grades introduced in our agent based model we were able to split transmission at school into three levels, so that a child can transmit to classmates, to students belonging to the same grade but different class, or to students from other grades, with specific transmission rates as detailed in [9]. In the absence of specific information we assumed transmission rates for class, grade and school to be independent of school level (i.e. primary or secondary education). Since pre-primary schools were not split into classes and grades and their average size was smaller than those of the other educational levels, the corresponding transmission rate was assumed equal to that in class for the other school levels. The model has therefore six transmission rates as a whole.
As regards infectivity over time, following state-of-the-art estimates on generation time distributions for the 2009 H1N1 influenza pandemic [9], we assumed infectiousness profiles over time to depend on the age of the case and on the setting where interaction occurs. We categorized the population into "students" and "adults" and assumed that interaction would occur in household, at school/workplace and in the general community. We assumed infectivity profiles to be discretized Weibull distributions with an offset of 1 (which is equivalent to having a latency period of 1 day), with parameters depending on whether the case was a student interacting in household, a student interacting outside the household (i.e. at school or in the community), or an adult. At any time , any infectious individual can infect others in every setting he belongs to, with probability ( ) = 1 − exp (− Δt), where ( ) is the force of infection exerted at time by individual in that specific setting. For instance, the force of infection of infectious individual on another member of the same household is is the infectivity of the case in the household, which is a function of the time of infection and of the age of the case, is the transmission rate in the household and is the household size. Analogous expressions describe the force of infection in schools, workplaces and general community. We assumed that, during school closure, household transmission remained unchanged, but community transmission increased of 25% [4].

Natural history of influenza and model calibration
Given our choice to consider infectivity profiles to depend on age (student/adult) and transmission setting, we assumed the value of the mean generation time to be 3.7 days (standard deviation 3.1 days) for students in household; 1.1 days (standard deviation 0.4 days) for students at school and in the community; 2.3 days (standard deviation 2.9 days) for adults [9].
The transmission model was parameterized in such a way that the overall proportions of cases in the different settings, assuming no differential susceptibility between adults and children, were consistent with empirical estimates, namely: 30% in households; 18% in schools; 19% in workplaces; 33% in the general community [4,[10][11][12], while assuming age-dependent susceptibility to infection led to 27% of transmission in households, 40% in schools, 8% in workplaces and 25% in the general community, similarly to what has been found in [2]. In addition, to match findings of [9], the proportion of transmission in class was set to be the same as the proportion of transmission in the entire school, and around 1.35 times the proportion of transmission in the grade. Different transmissibility scenarios were obtained by multiplying all transmission rates for suitable constant values.
Choosing generation time distributions and transmission rates determines the evolution of the epidemic.
Following the argument proposed in [13] the basic reproduction number 0 is generally defined as the average number of secondary infections caused by a typical primary infection in a fully susceptible population [14]. However, the definition of a typical infectious individual for a model with highly heterogeneous mixing such as the one used in this work is not straightforward. Therefore, in the parameterization procedure we considered in place of 0 the quantity , which is defined as the average number of secondary infections generated by a randomly selected individual in a fully susceptible population [4]. Ferguson et al. [4] showed that for US and UK populations the relation 0 = + 0.2 holds. In our simulations yielding = 1.5 we estimate ad epidemic growth rate of 0.221 days -1 , in excellent agreement with the estimate of 0.219 days -1 in [4] for simulations having = 1.5 and 0 = 1.7. This allows presenting all results in terms of 0 rather than , making it easier to compare them with the ones reported in the literature.
For classic homogeneous mixing SIR models, it is well known [4,15] that 0 = 1 + , where is the overall generation time and is the exponential growth rate of the epidemic. Since we are considering setting and age-specific generation time distributions, we do not know the value of the overall generation time for our model; moreover, we want to test whether this approximation still holds given the high heterogeneity of our model and, if so, for which value of the generation time (named ). By minimizing the quantity where is the growth rate of the simulated epidemic having 0 = , we obtain = 3.2 days, in perfect agreement with the findings reported in [4] and we can conclude that the relation between 0 and is a good approximation in the explored range 0 = 1.1, … ,2.3. This justifies the choice, already used for instance in [2,8], of multiplying transmission rates for suitable scaling factors in order to obtain different transmissibility scenarios.

Basic reproductive number
The transmission model was calibrated assuming a basic reproductive number R0=1.5, to be consistent with estimates for the 2009 H1N1 influenza pandemic presented in [16]. R0 was computed from the intrinsic growth rate during the initial phase of the simulated epidemic without assuming any kind of intervention, initialized with 100 infectious individuals randomly chosen in the population so that stochastic effects were partially eliminated. As a sensitivity analysis we also explored scenarios with R0=1.3, 1.7.

Probability of being symptomatic
In the reference scenario we considered the proportion of influenza cases showing symptoms to be 30%, as described in [17]. 50% probability of being symptomatic was also explored as sensitivity analysis.

Differential susceptibility by age
In the reference scenario we assumed children to be twice more susceptible to infection than adults, as estimated for the 2009 A/H1N1 influenza pandemic as well as for other influenza epidemics [9,18,19]; in the main text we also explored the scenarios where children are either as susceptible or four times more susceptible than adults.

Time of simulation
We set the day when the first case in the population occurred as the initial time of the simulation, and let the simulation run for 600 days. Closure interventions can be implemented as soon as the cumulative number of cases experiencing symptoms in the UK as a whole reaches 0.25% of the total population (corresponding to ~140,000 individuals), and during the following 180 days.

Number of realizations
In order to ensure stability of findings, all presented results were obtained by averaging over 50 stochastic realizations of the same experiment. Around 26,000 simulations were carried out as a whole to produce the full analysis. Variability among simulations was larger for the epidemic timing, due to the stochastic nature of the initial phases, while final infection attack rate, depending on the course of the entire epidemic, was very stable.

Start of interventions
The choice, made in the main text, to begin the application of interventions as soon as the cumulative number of cases experiencing symptoms in the UK as a whole reaches 0.25% or higher of the total population is arbitrary. To perform a sensitivity analysis on this threshold, we considered the framework corresponding to the best attack rate reduction provided to lose 4 weeks, and assumed to allow closures after the first case in the population (which represents an extreme situation) or as soon as cumulative symptomatic cases were 0.5% or 1% of the total population. We obtained that, for all closure types, all thresholds on the initial cumulative incidence yield attack rate reduction comparable to the baseline 0.25% threshold, while starting closures after the very first case would be too early, causing too little impact on attack rate (Fig S1). This means that a cumulative case incidence of 0.25% or higher may be a good trigger for starting closures while avoiding false alarms that may not be unusual in the very first phases of the epidemic.

Transient and permanent effects
All results presented in the main text are evaluated right after the end of the period during which application of closure policies is possible: we may refer to this as the "transient" behavior of the epidemic. Of course, in principle the epidemic may not be over by that time and a resurgence of cases would still be possible, so we analyze the main epidemiological indicators (i.e. attack rate reduction, peak incidence reduction, peak delay) of the best strategies also at the end of the simulation. We assumed the final time of each simulation to be 600 days after the first national case, without considering holidays typically planned during the school year in order to emphasize the effect of unscheduled closures. We infer (Fig S2) that the only situation where some significant difference between transient and final behavior can be observed is attack rate reduction for 8 weeks lost: in that case best reactive, gradual and county closures yield a 40-45% attack rate reduction with respect to the baseline scenario at the end of closure implementation, dropping to about 30% at the end of the simulation. The only other marginal difference between transient and final attack rate reduction can be observed for county closure at 4 weeks lost. This means that all best implementations (for whatever indicator and closure type), applicable in a period lasting 180 days, are almost always able to halt the epidemic, with the only relevant exception of 8 weeks lost, where a second wave of cases at the end of the closure implementation period can be observed (see Fig S3).  Fig S3. Daily incidence. National daily cases per 10,000 individuals for the epidemic without any closure intervention (grey) and epidemic resulting from optimal closure strategies in terms of transient attack rate reduction (red: reactive closure; green: gradual closure; blue: county closure). Solid lines represent mean values from the total simulation runs; shaded regions represent the corresponding 95% confidence intervals. Vertical lines correspond to the starting and ending dates of the period during which closures can be implemented, with their 95% confidence intervals represented on top. a, c, e. Four weeks lost per student. b, d, f. Eight weeks lost per student.