Impacts of K-12 school reopening on the COVID-19 epidemic in Indiana, USA

Graphical abstract


Introduction
The United States has been the country most severely impacted by the COVID-19 pandemic in terms of total reported cases and deaths, with over 28 million reported cases and more than 500 thousand deaths by March, 2021(The New York Times, 2020. This severity led to social interventions on an unprecedented scale, including restrictions on mass gatherings, bans on non-essential travel, and school closures Flaxman et al., 2020;Coronavirus Government Response Tracker, 2021;Badr et al., 2020). While such restrictions were initially successful in reducing transmission, the subsequent relaxation of restrictions on mass gatherings and movement were followed by large increases in notified cases and deaths (The New York Times, 2020;Flaxman et al., 2020;Auger et al., 2020;Courtemanche et al., 2020). By the time the 2020-2021 school year began in August, transmission was at its highest point in the epidemic yet in some parts of the US. In Indiana, for example, the maximum number of daily cases was around 1, 200 by then, which was higher than the previous maximum of fewer than 800 in late April (ISDH, 2021).
This context of intense community transmission raised numerous questions about how schools should approach reopening for the start of the school year in August (Dibner et al., 2020; The White House, 2021; Indiana State Department of Health and Indiana Department of Education, Indiana Family and Social Services Administration, 2020). During influenza epidemics, school closures have been estimated to reduce transmission community-wide (Jackson et al., 2013;Kawano and Kakehashi, 2015;Cauchemez et al., 2009). In general, schools are seen as key drivers of the transmission of respiratory pathogens due to close contact among children at school (Hens et al., 2009;Cauchemez et al., 2008;Jackson et al., 2014). However, several factors complicated the effect of school reopenings on SARS-CoV-2 transmission. In particular, children and adolescents appear less susceptible to infection and are much less likely to experience severe outcomes following infection (Vermund and Pitzer, 2020;Davies et al., 2020;Zhang et al., 2020;Wu et al., 2020;Verity et al., 2020;Willem et al., 2020). It is also still unclear what their contribution to transmission is, but several studies suggest they can play an important role (Vermund and Pitzer, 2020;Szablewski, 2020;Park et al., 2020;Fateh-Moghadam et al., 2020). There are additional economic and social factors to consider, too, such as the economic costs of school closures for families that must then stay home from work, and the nutritional benefits of school reopening for children who rely on free and subsidized school meals (Psacharopoulos et al., 2020;Dooley et al., 2020;Van Lancker and Parolin, 2020).
Our objective in this study was to explore how different conditions for school reopening during the fall semester of 2020 could have impacted the statewide burden of COVID-19 in Indiana. Specifically, we focused on the effects of school operating capacity and adherence to wearing face masks in schools. This focus was motivated by the fact that Indiana and other US states reopened their schools for in-person instruction in August with only minimal interventions of requiring face masks and physical distancing in schools, despite uncertainty about the proportion of students who would elect to attend in person and the degree to which they would adhere to face-mask and physical-distancing requirements. We approached this question with an agent-based model originally developed for pandemic influenza (Grefenstette et al., 2013), which we tailored to SARS-CoV-2 (Davies et al., 2020;Wu et al., 2020;Bi et al., 2020;Lauer et al., 2020;Ganyani et al., 2020;Chu et al., 2020) and applied to a geographically and demographically realistic synthetic population representing Indiana. In addition to presenting outcomes across a range of hypothetical scenarios, we calibrated our model to data from the fall to assess the plausibility of K-12 school reopening as a driver of the observed resurgence of SARS-CoV-2 in Indiana during fall 2020.

Approach
Our approach to modeling SARS-CoV-2 transmission was based on the Framework for Reconstructing Epidemic Dynamics (FRED) (Grefenstette et al., 2013), an agent-based model that offers the ability to explore the impacts of complex, non-pharmaceutical interventions in a natural way through modifications to individual behaviors. Using this model, we simulated the spread of SARS-CoV-2 in Indiana using a synthetic population with demographic and geographic characteristics of the state's real population, including age, household composition, household location, and occupation (Wheaton, 2012). We analyzed the impact of school reopening from August 24 (first day of classes in Marion County, the most populous) to December 31, 2020 in the population of Indiana as a whole, as well as in students, teachers, and their cohabitants. We quantified impact as the difference in the number of COVID-19 infections, symptomatic infections, and deaths between each scenario and a baseline scenario, the details of which differed according to which comparisons were of interest.

Agent-based model
Our model was based on a synthetic population of the entire state of Indiana (Wheaton, 2012), which included 1.3 million K-12 students, 1.7 million people living with students, and 6.3 million people in total (Table S10). Each of these agents visits a set of places defined by their activity space, which can include houses, schools, workplaces, long-term care facilities, and various neighborhood locations. Transmission can occur when an infected agent visits the same location as a susceptible agent on the same day, with numbers of contacts per agent specific to each location type. For example, school contacts depend not on the size of the school but on the age of the student and their assigned school grade, given that students have a higher number of contacts with students in their classroom than with those in other classrooms. Every day of the week, students and teachers visit their school, and students are assigned to classrooms based on their age. Given that schools are closed during the weekends, community contact is increased by 50 % for students and teachers on weekends, unless they are sheltering in place (Grefenstette et al., 2013). Community contacts are modeled as contacts with other agents in the same neighborhood. For both schools and other locations, we adopted contact rates for each location type that were previously calibrated to influenza attack rates specific to each location type (Grefenstette et al., 2013;Mossong et al., 2008).
Once infected, each agent had a latent period (mean = 3.35 days, standard deviation = 1.16 d) and an infectious period (mean = 3.7 d, s.d. = 1.2 d) drawn from distributions calibrated so that the average generation interval distribution matched estimates from Singapore (mean = 5.20 d, s.d. = 1.72 d) (Ganyani et al., 2020). The absolute risk of transmission depended on the number and location of an infected agent's contacts and a parameter that controls SARS-CoV-2 transmissibility upon contact, which we calibrated. A proportion of the infections were asymptomatic, with the probability of symptoms increasing with age (Davies et al., 2020;Mizumoto et al., 2020;The Novel Coronavirus Pneumonia Emergency Response Epidemiology Team, 2020;Emery et al., 2020) (Table S9, Fig. S4). We assumed that these infections were as infectious as symptomatic infections and had identical incubation and infectious period distributions (Fateh-Moghadam et al., 2020;Yonker et al., 2020;Heald-Sargent et al., 2020;Hu et al., 2020). Furthermore, we assumed that children were less susceptible to infection than adults, which we modeled with a logistic function calibrated to model-based estimates of this relationship by (Davies et al., 2020). We assumed that severity of disease increased with age, consistent with statistical analyses described elsewhere Bi et al., 2020;Lauer et al., 2020;The Novel Coronavirus Pneumonia Emergency Response Epidemiology Team, 2020).
Agent behavior in FRED has the potential to change over the course of an epidemic. Following the onset of symptoms, infected agents selfisolate at home according to a fixed daily probability, whereas others continue their daily activities (Perkins et al., 2016;COVID T.C. et al., 2020). This probability was chosen so that, on average, 68 % of agents will self-isolate at some point during their symptoms. This figure of 68 % was based on the proportion of symptomatic infections among healthcare workers in the USA that developed fever during the course of their infection (COVID T.C. et al., 2020), and assuming that those with fever are likely to self-isolate. Agents can also engage in a variety of non-pharmaceutical interventions, including school closure, sheltering in place, and a combination of mask-wearing and physical distancing. School closures occur on specific dates across the state (State of Indiana, 2021), resulting in students limiting their activity space to household and neighborhood locations on those days. Within households, agents interacted with their cohabitants on a daily basis. We assumed that agents did not wear face masks inside their homes, nor did they isolate from their household members if infected. To capture temporal changes in overall mobility and community contact over the course of the epidemic, we modeled a time-varying probability of sheltering in place. On days when an agent sheltered in place, they reduced their activity spaces to their home only, whereas other agents continued with their normal routines.
We modeled protection from face masks by reducing the probability of transmission when an infected agent wore a face mask. Our default assumption followed a median estimate of an adjusted odds ratio of 0.3 against SARS-CoV in non-healthcare settings (Chu et al., 2020;Wu et al., 2004). A meta-analysis by Chu et al. (2020) included studies in healthcare settings in addition to those from non-healthcare settings, but we included only estimates that referred to the latter. The proportion of agents wearing face masks in workplace and community settings changed over time, and we assumed that agents did not wear face masks inside their households. In addition, students and teachers wore face masks according to probabilities specified as part of school reopening scenarios in our analysis. Further details about the model are available in the Supplementary Text.

Model calibration and validation
Two of the time-varying drivers of transmission in our model were informed by time-varying data inputs. First, we informed the daily probability of sheltering in place with mobility reports from Google (COVID-19 Community Mobility Report, 2021). In doing so, this sheltering-in-place probability in our model accounts for both the effects of shelter-in-place orders and some people deciding to continue staying at home after those orders are lifted (Cronin and Evans, 2020). Second, we informed the daily proportion of agents wearing face masks in workplace and community settings with Google Trends data for Indiana using the terms "face mask" and "social distancing" (Google Trends, 2021). To inform the magnitude of the proportion of people wearing face masks in workplace and community settings, we used survey data (New York Times and Dynata, 2021) on face-mask usage from a single point in time.
The values of nine model parameters (listed in Table S9) were informed by calibrating the model to three time-varying epidemiological data streams corresponding to the state of Indiana-daily deaths, daily hospitalizations, and daily test positivity at the state level through August 10-and to the age distribution of cumulative deaths through July 13. We obtained daily incidence of reported cases and deaths from the New York Times COVID-19 database (The New York Times, 2020). Daily hospitalizations and the age distribution of cumulative deaths were obtained from the Indiana COVID-19 dashboard (ISDH, 2021). Daily numbers of tests performed in the state were available from The Covid Project (The COVID Tracking Project, 2020). To calibrate the model to these data, we first used a Sobol design sampling algorithm (King et al., 2020(King et al., , 2016 to draw 6,000 combinations of the nine calibrated parameters. We then calculated the likelihood of each parameter combination given the data, and resampled the parameter combinations proportional to their likelihoods to obtain an approximation of the posterior distribution of parameter values. Additional details about the calibration procedure are described in the Supplementary Material. To validate the model, we compared its predictions to data withheld from model calibration. Specifically, we compared the calibrated model's predictions of the infection attack rate, both overall and by age, to results from two statewide serological surveys undertaken in late April and early June (Menachemi et al., 2020;Zeek, 2020). We assessed the model's success in this validation exercise by visual comparison of the model's predictions with the data, with a focus on overlap between the 95 % prediction intervals of the model and the 95 % credible intervals associated with the empirically-derived estimates from the serological surveys.

Model outputs
The main outputs from our model were the numbers of infections, symptomatic cases, hospitalizations, and deaths at the state level. For different comparisons, these outputs were examined either on a daily basis, cumulatively between August 24 and December 31, or stratified by place of infection (school, home, other) or affiliation with schools (student, teacher, none). Another output that we examined was the daily reproduction number, R(t), which was calculated in the model for each infected agent and averaged across the modeled population. We defined R(t) as the average number of secondary infections caused by any agent who was infected on day t -i.e., the case reproduction number a la Fraser (Fraser, 2007). Finally, to account for uncertainty in these outcomes due to uncertainty in the calibrated parameters, we sampled parameter sets from our approximation of the posterior distribution of parameter values.

Model scenarios 2.5.1. Effects of conditions in schools
To explore how alternative conditions for school reopening could have impacted the statewide burden of COVID-19 in Indiana, we performed simulations that spanned a range of assumptions about school operating capacity and adherence of students and teachers to face masks while in school. We chose to focus on these parameters given that they were two of the major unknowns as the state proceeded with its plans for in-person instruction beginning in August 2020. Given that students were offered the option of either in-person or remote instruction, we evaluated scenarios in which school operating capacity was either 50 %, 75 %, or 100 %. Specifically, this parameter represents the daily probability that a student would go to school, such that all students could go to school at some point during the simulation but the average number of students in attendance on a given day is determined by the school operating capacity parameter. For each of these scenarios about school operating capacity, we also considered scenarios in which the adherence of students and teachers to face masks was either 50 %, 75 %, or 100 %. In addition, we considered a scenario in which schools reopened normally (100 % capacity, 0 % face-mask adherence) and a scenario in which schools operated remotely (0 % capacity, face-mask adherence irrelevant since no one in school).

Sensitivity analyses
To explore the sensitivity of our results to model uncertainties, we performed three sets of sensitivity analysis. First, for each of the nine scenarios comprising our primary analysis, we analyzed the sensitivity of cumulative infections and the proportion of infections acquired in schools to each of the nine calibrated parameters by calculating partial rank correlation coefficients (Nunes MS with contributions from T et al., 2020). Second, we considered the sensitivity of our results to values of two assumed parameters not included in the calibration: protection afforded by face masks and the probability of isolation given symptoms. This included a total of four scenarios exploring lower and higher values of each of those parameters. Third, we considered alternative scenarios about select model assumptions that we regarded as potentially important unknowns about transmission of SARS-CoV-2 by children. These included scenarios in which asymptomatic infections (which are more likely to occur among children) are half as infectious as symptomatic infections, a scenario in which children aged 0-10 years have lower susceptibility (0.1), and a scenario in which individuals of all ages are equally susceptible to SARS-CoV-2 infection. For each of the alternative scenarios in the second and third sets of sensitivity analyses, we re-calibrated the model under that scenario and simulated it forward for the fall semester under the nine primary scenarios about school operating capacity and adherence to face-masks in schools. These latter two sets of sensitivity analyses all focused on cumulative infections statewide between August 24 and December 31, 2020.

Retrospective analysis
To understand what conditions for school reopening resulted in model behavior consistent with data from fall 2020, we calibrated parameters for school operating capacity and face-mask adherence in schools to data from this time period. To do so, we held all other calibrated parameters at the values we estimated based on calibration to data from January 1 to August 24, 2020. Whereas we used Google mobility reports to drive a time-varying probability of sheltering in place during that initial period, we opted to assume fixed levels of mobility from August 24 through December 31. The main reason for this choice was that the Google data we used for the initial period showed a decrease in mobility during the fall, despite an increase in incidence, making it difficult to explain the epidemiological data under that assumption. To overcome that problem and to make this analysis more directly comparable to our other analyses, we held the probability of sheltering in place fixed at either of two levels of mobility during the fall: 1) mobility remained at summer levels, or 2) mobility increased to prepandemic levels. Under these assumptions, we calibrated the two focal parameters about conditions in schools to data from the fall using the same calibration procedure as we used during the initial period.

Model calibration and validation
Our model was generally consistent with the data to which it was calibrated, capturing trends over time in daily deaths, hospitalizations, and test positivity at the state level ( Fig. 1A-C), as well as greater proportions of deaths among older age groups (Fig. 1D). Some trade-offs in the model's ability to recreate different data types were apparent, such as a recent increase in hospitalizations that the model failed to capture (Fig. 1C), likely due to the predominance of data on deaths in the likelihood. Similarly, the model underestimated the proportion of deaths in people older than 80 (Fig. 1B), indicating a possible underestimation of the contacts in this age group in the overall community or in long-term care facilities. Even so, the model's predictions reproduced the range of variability in the data, as assessed by the coverage probabilities of its 95 % posterior predictive intervals (daily deaths: 0.85; daily hospitalizations: 0.93; daily test positivity: 0.95; cumulative deaths by age: 1.0). The model was also consistent with data withheld from fitting. Across all ages, the model's 95 % posterior predictive intervals of the cumulative proportion infected through late April (median: 0.017; 95 % CrI: 0.0045− 0.051) and early June (median: 0.022; 95 % CrI: 0.0058− 0.069) spanned estimates from two state-wide serological surveys (Menachemi et al., 2020) (Fig. 2A). Our model's predictions also overlapped with age-stratified estimates from those surveys (Fig. 2B), although it underpredicted infections among individuals aged 40-60 years.
Calibration of the parameter that scaled the magnitude of SARS-CoV- Fig. 1. Model calibration to statewide data: A) daily incidence of death; B) proportion of deaths through July 13 in decadal age bins; C) daily incidence of hospitalization; and D) daily proportion of tests administered that are positive for SARS-CoV-2. In all panels, blue diamonds represent data. In A, C, and D, the gray line is the median, the dark shaded region the 50 % posterior predictive interval, and the light shaded region the 95 % posterior predictive interval.
2 importations (Perkins et al., 2020;MIDAS Network, 2021) (Fig. 1A). This estimate is within the range of other estimates for R(t) that include the state of Indiana (Fernández-Villaverde and Jones, 2020; Mohler et al., 2021). Driven by a calibrated estimate that the proportion of people sheltering in place rose in early March and peaked at a median of 32.1 % (95 % CrI: 28.8-66.9 %) on April 7 (Fig. S2A), our estimates of R(t) dropped to a low of 0.57 (95 % CrI: 0.42− 0.71) on April 7 and remained below 1 thereafter (Fig. 1A). Also impacting our estimates of R(t) was the increasing use of face masks in the community, which we estimated at 53.4 % (95 % CrI: 46.1-54.0 %) as of July 19 (Fig. S2B). Note that this estimated distribution of community face-mask adherence does not differ between scenarios and is not the same as the level of adherence in schools, which we imposed at different levels depending on the scenario.

Effects on statewide burden
Under a scenario in which schools reopened at full capacity and without any use of face masks, our model projected that R(t) across the state as a whole would have increased to 1.72 (95 % CrI: 1.43-2.17) by mid-September (Fig. 3A). Given our assumption that levels of sheltering in place and face-mask adherence in the community remained constant during the fall (Fig. S2B), this increase in transmission was driven by infections arising in schools (Fig. 3B). As a result, new infections statewide would have risen to levels in the fall far exceeding those from the spring (Fig. 3C). Under this scenario, our model projected a total of 2.57 million (95 % CrI: 2.36-2.88 million) infections (Fig. 4A)  Under a scenario in which schools went to remote instruction and all children remained at home, our model projected that R(t) would have remained near levels from August for the remainder of 2020 (Fig. S3A). This was a result of our assumption that sheltering in place and facemask adherence in the community would have remained at their estimated levels as of August 13. Under this scenario, transmission would have continued through contacts at workplaces, within homes, and elsewhere in the community (Fig. S3B Less extreme scenarios about school operating capacity and facemask adherence in schools also resulted in a wide range of variation in the projected statewide burden of COVID-19 during fall 2020. Under a scenario in which schools operated at 50 % capacity and achieved 100 % face-mask adherence, the cumulative numbers of infections and deaths that our model projected were similar to projections under a scenario in which schools operated remotely (Fig. 4, Tables S1 & S5). In general, cumulative infections and deaths statewide in fall 2020 were more sensitive to school operating capacity than to face-mask adherence in schools, with the worst outcomes projected to occur under a scenario with 100 % school operating capacity and 50 % face-mask adherence. Under this scenario, cumulative infections statewide were projected to have been 42.8 (95 % CrI: 41.3-44.3) times greater than if schools had operated remotely, and cumulative deaths statewide were projected to have been 9.2 (95 % CrI: 8.9-9.5) times greater (Table S1).

Effects on risk for individuals affiliated with schools
Relative to a scenario with remote instruction, risk of infection and symptomatic infection was greatest for students (Figs. 5 & S5, left column), with a hundred-fold or greater increase in the risk of infection under a scenario with 100 % school operating capacity and 50-75 % face-mask adherence in schools (Tables S2 & S6). The risk of symptomatic disease in school-aged children was two-fold lower for children under 10 years of age (Fig. S6). Compared to students, the risk of infection was slightly lower for teachers, and much lower for students' families. Due to their older ages, however, teachers and families experienced a much higher risk of death than students (Figs. 5 & S5, center & right columns). The highest risk of death was for teachers under scenarios with 100 % school operating capacity (Figs. 5 & S5, center column). Compared to a scenario with remote instruction, the relative risk of death for teachers under these scenarios ranged from a 41-fold increase when face-mask adherence was 100 % to a 166-fold increase when face-mask adherence was 50 % (Table S3, S7). Under scenarios with 75 % school operating capacity, those same relative risks dropped to a four-fold increase when face-mask adherence was 100 % and a 22fold increase when face-mask adherence was 50 %. This again illustrates the overall greater effect of school operating capacity than face-mask adherence in schools.

Sensitivity analyses
Our sensitivity analysis of the model's nine calibrated parameters quantified the partial rank correlation coefficient (PRCC) of each of two model outputs: cumulative infections statewide from August 24 to December 31, and the proportion of infections acquired in schools. In general, these outputs were most sensitive to parameters controlling the age-susceptibility relationship (Figs. S7 & S8). Under some scenarios, the minimum susceptibility parameter, which applied to young children, had a PRCC as high as 0.6. The transmissibility parameter also had a PRCC that high, but only in some scenarios. For example, in scenarios with lower school operating capacity, the two parameters most relevant to community transmission-transmissibility and face-mask adherence in community settings-had a greater influence on cumulative infections statewide (Fig. S7), given that the contribution of schools to transmission was diminished in those scenarios.
Our sensitivity analysis of parameters for protection afforded by face masks and the probability of isolation given symptoms focused on cumulative infections statewide in fall 2020. For each alternative value of these parameters, we re-calibrated the model and generated projections under our nine primary scenarios about school operating capacity and face-mask adherence in schools. In general, the relative effects of differences in face-mask adherence in schools and school operating capacity were similar under alternative values of protection afforded by face masks (Table S11). However, relative to a baseline with schools operating remotely, the magnitude of the proportional increase in cumulative infections was sensitive to the level of protection afforded by face masks. For example, under a scenario with 75 % face-mask adherence and 75 % school operating capacity, the increase in cumulative infections ranged from 23.5-fold to 1.6-fold across the range of values of protection afforded by face masks that we explored (adjusted odds ratio = 0.12− 0.73). Proportional increases in cumulative infections were generally insensitive to the probability of isolation given symptoms (Table S12). For example, under a scenario with 75 % facemask adherence and 75 % school operating capacity, the increase in cumulative infections ranged from 4.9-fold to 6.9-fold across the range of values of isolation probability that we explored (0.5− 0.9).
Our sensitivity analysis of assumptions related to the role of children in transmission also focused on cumulative infections statewide in fall 2020. Relative to a baseline with schools operating remotely, the proportional increase in cumulative infections was very similar to our default assumptions under a scenario with lower susceptibility among children aged 0-10 years and a scenario with lower infectiousness of asymptomatic infections (Table S13). This was the case for all scenarios about face-mask adherence in schools and school operating capacity, except for the most extreme case in which face-mask adherence in schools was 50 % and school operating capacity was 100 %. In that case, our default assumptions resulted in a 42.9-fold increase in cumulative infections, whereas the two alternative scenarios resulted in a 26-fold increase. Under a scenario with equal susceptibility for all ages, proportional increases in cumulative infections were much higher than under our default assumptions (Table S13). Because susceptibility in children was higher under this scenario, transmission statewide was much higher when school reopened, especially in scenarios with higher school operating capacity (Fig. S9). By the same token, prevalence dropped to levels over the summer when school was not in session that were much lower than the data suggest (Fig. S10), which raises doubts about the plausibility of this scenario.

Retrospective analysis
The model successfully reproduced statewide data from fall 2020 under relatively high values of school operating capacity and intermediate values of face-mask adherence in schools (Fig. 6). Under an assumption that the probability of sheltering in place in the state as a whole was fixed at levels from summer 2020 (Fig. 6, gray), the model calibration resulted in a median school operating capacity of 0.88 (95 % CrI: 0.8− 0.95) and a median face-mask adherence in schools of 0.6 (95 % CrI: 0.4− 0.7). Under an assumption that the probability of sheltering in place was fixed at pre-pandemic levels (Fig. 6, red), the model calibration resulted in a median school operating capacity of 0.85 (95 % CrI: 0.77− 0.92) and a median face-mask adherence in schools of 0.61 (95 % CrI: 0.41− 0.75). Overall, there was a wide range of values of face-mask adherence in schools that were consistent with the data, and the range of values of school operating capacity consistent with the data depended somewhat on face-mask adherence in schools (Fig. 6, bottom left).

Discussion
Our model provides a detailed, demographically realistic representation of SARS-CoV-2 transmission in Indiana that is consistent both with data to which it was calibrated and to data that was withheld from calibration. In contrast to models that rely on assumptions about intervention impacts or estimate them statistically Esposito and Principi, 2020), our model makes predictions about intervention impacts based on first-principles assumptions about individual-level behavior and contact patterns. Consistent with results from other analyses Zhang et al., 2020;Esposito and Principi, 2020), the inputs and assumptions in our model led to a prediction that schools made a considerable contribution to SARS-CoV-2 transmission in February and early March, prior to large-scale changes in behavior. Extending that, a primary result of our analysis is that K-12 school reopening was capable of making a considerable contribution to SARS-CoV-2 transmission during fall of 2020, with the degree of that contribution dependent on conditions in schools.
The burden of COVID-19 associated with in-person school operation was predicted by our model to fall unevenly across the state's population. In scenarios of school operating capacity of 100 % with 50 % face-mask adherence in schools, our model predicted that hundreds of thousands of children could have been infected during the fall semester, with very few of those resulting in deaths. In contrast, our results show that hundreds of deaths in teachers and school-affiliated families could have occurred. Our model indicates that the burden of COVID-19 in , and deaths (bottom row) per 1,000 people. These outcomes are presented separately for students (left column), teachers (middle column), and school-affiliated families (right column). Scenarios are defined by school operating capacity (x-axis) and face-mask adherence in schools (shading). Orange lines represent projections under a scenario of school reopening at full capacity without masks (solid: median; dotted: 95 % posterior predictive interval). Blue lines represent a scenario where schools operate remotely. Error bars indicate inter-quartile ranges. schools, teachers, and school-affiliated families across the state could have been reduced by operating at reduced capacity and achieving high face-mask adherence in schools. Under the relatively optimistic scenario of 50 % school operating capacity and 100 % face-mask adherence in schools, our model predicted that infections and deaths statewide would have been only 22 % greater than under a scenario with fully remote instruction. In contrast, our model results suggest that if schools would have operated at full capacity, infections and deaths statewide could have been one to two orders of magnitude greater than the scenario with fully remote instruction, especially with poor face-mask adherence in schools. When we extended our model calibration to account for data from fall 2020, we found that school operating capacity of 80-95 % and face-mask adherence in schools of 40-70 % resulted in model predictions most consistent with the observed data. For reference, data from the National COVID-19 School Response Dashboard (COVID-19 Resource Hub: STATS Indiana, 2021) indicate that school operating capacity in Indiana was 77-83 % during the fall. Although conditions elsewhere in the community likely played a role in statewide trends during the fall and were not accounted for fully by our model, these results demonstrate that transmission associated with K-12 school reopening was capable of driving the statewide resurgence of COVID-19 observed in Indiana in fall 2020.
The impacts associated with reduced school operating capacity result from reductions in both the number of contacts within the school and the probability that an infected student would be in attendance in the first place, similar to the logic behind why smaller gatherings are associated with reduced risk of transmission (Flaxman et al., 2020;Che Mat et al., 2020;Jang et al., 2020). The magnitude of our results was most sensitive Fig. 6. Retrospective analysis of the model calibrated to statewide data from fall 2020. Under two alternative scenarios about the daily probability of sheltering in place (summer vs. pre-pandemic mobility level in gray and red, respectively), we calibrated the parameters for school operating capacity and face-mask adherence in schools to data from August 24 through December 31, 2020 (blue diamonds). The calibrated model's correspondence to daily incidence of death statewide is shown in the top two panels, and values of the calibrated parameters are shown in the bottom panels.
to the degree of protection afforded by face masks, which remains uncertain in school and other community settings for SARS-CoV-2 (Chu et al., 2020). That uncertainty can be reduced as more studies are conducted. Recently, some studies have shown that face masks offer significant protection in community settings similar to what we assumed. For example, Payne et al. followed 382 U.S. Navy service members who reported wearing face masks and found a reduced risk of 70 % in those who did (Payne et al., 2020). Similar values were found in studies of 124 households in China  and 839 close contacts of 211 index cases in Thailand (Doung-Ngern et al., 2020).
Although the scenarios we considered resulted in projected impacts spanning nearly the full range between fully remote instruction and fully in-person instruction with no face masks, they are a simplification of the complexities of how schools likely operated in fall 2020. Scenarios that we did not explore include different groups of students attending in person or remotely (Keeling et al., 2020), varying degrees of modularization within schools (Head et al., 2020), and the implementation of testing-based control strategies in schools (Panovska-Griffiths et al., 2020). In the event that infectiousness is lower for asymptomatic infections, the impact of school reopening on lower grades could have been lower than our results suggest. A related simplification of our statewide analysis is that the state, in reality, consists of a patchwork of policies across districts. In light of this complexity that our model does not capture, our results should be interpreted with caution in relation to specific counties or school districts below the state level. Across all scenarios though, our results illustrate the importance of reduced school operating capacity and maximal face-mask adherence in schools, as do other modeling studies (Keeling et al., 2020;Head et al., 2020;Panovska-Griffiths et al., 2020;Lee et al., 2020;Landeros et al., 2020).
A critical assumption of our analysis is that children are capable of being infected with SARS-CoV-2 and transmitting it to others at meaningful levels. Although the burden of severe disease skews strongly towards older ages (ISDH, 2021;Verity et al., 2020;COVID-19 Hospitalizations, 2021), there are other lines of evidence that support our assumption. These include a contact-tracing study that found no distinguishable difference between infectivity of children and adults (Fateh-Moghadam et al., 2020), several studies that found no distinguishable difference in viral load between children and adults (Yonker et al., 2020;Heald-Sargent et al., 2020;Lavezzo et al., 2020;Jones et al., 2020), a study that observed a greater secondary attack rate among children in homes (Fateh-Moghadam et al., 2020), and a modeling study that found no evidence that children were less infectious (Dattner et al., 2020). More direct evidence comes from COVID-19 outbreaks that have been observed in schools, such as one in a high school in Israel in which 13.2 % of students and 16.6 % of staff were infected in just 10 days (Stein-Zamir et al., 2020). Even more pertinent, since schools reopened in September 2020, 31,658 COVID-19 cases have been reported in students and 13,240 cases have been reported in teachers and school staff across the state, as of April 2021 (ISDH, 2021).
There is now a growing body of evidence that school closures contributed to mitigating the first wave of the epidemic (Auger et al., 2020;Panovska-Griffiths et al., 2020) and, as we have shown, may have contributed to the resurgence of SARS-CoV-2 during fall 2020. Our study adds to this evidence, and suggests an even greater impact of school reopening than several other studies (Keeling et al., 2020;Head et al., 2020;Panovska-Griffiths et al., 2020;Landeros et al., 2020;Abdollahi et al., 2020). This is due in part to our assumption that asymptomatic and symptomatic infections contribute similarly to transmission (Fateh-Moghadam et al., 2020;Yonker et al., 2020;Heald-Sargent et al., 2020;Lavezzo et al., 2020;Jones et al., 2020), and in part to our model's ability to capture chains of transmission within schools and extending out into the community. Our study echoes several modeling studies in emphasizing the importance of reducing school operating capacity to impede transmission (Head et al., 2020;Panovska-Griffiths et al., 2020;Lee et al., 2020;Landeros et al., 2020;Abdollahi et al., 2020). As schools grapple with COVID-19 going forward, results such as these provide an important basis for motivating the adoption and sustainment of reduced school operating capacity and adherence to face-mask requirements in schools. As we demonstrated, these actions are highly consequential for those directly linked to schools and for the communities in which they are embedded.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.