Estimating the impact of school closures on the COVID-19 dynamics in 74 countries: A modelling analysis

Abstract

Background: School closures have been a prominent component of the global Coronavirus Disease 2019 (COVID-19) response. However, their effect on viral transmission, COVID-19 mortality and health care system pressure remains incompletely understood, as traditional observational studies fall short in assessing such population-level impacts.

Methods and findings: We used a mathematical model to simulate the COVID-19 epidemics of 74 countries, incorporating observed data from 2020 to 2022 and historical school closure timelines. We then simulated a counterfactual scenario, assuming that schools remained open throughout the study period. We compared the simulated epidemics in terms of Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) infections, deaths, and hospital occupancy pressure. We estimated that school closures achieved moderate to significant burden reductions in most settings over the period 2020 to 2022. They reduced peak hospital occupancy pressure in nearly all countries, with 72 out of 74 countries (97%) showing a positive median estimated effect, and median estimated effect ranging from reducing peak hospital occupancy pressure by 89% in Brazil to increasing it by 19% in Indonesia. The median estimated effect of school closures on COVID-19 deaths ranged from a 73% reduction in Thailand to a 7% increase in the United Kingdom. We estimated that school closures may have increased overall COVID-19 mortality (based on median estimates) in 9 countries (12%), including several European nations and Indonesia. This is attributed to changes in population-level immunity dynamics, leading to a concentration of the epidemic during the Delta variant period, alongside an upward shift in the age distribution of infections. While our estimates were associated with significant uncertainty, our sensitivity analyses exploring the impact of social mixing assumptions revealed robustness in our country-specific conclusions. The main study limitations include the fact that analyses were conducted at the national level, whereas school closure policies often varied by region. Furthermore, some regions, including Africa, were underrepresented due to insufficient data informing the model.

Conclusions: Our analysis revealed nuanced effects of school closures on COVID-19 dynamics, with reductions in COVID-19 impacts in most countries but negative epidemiological effects in a few others. We identified critical mechanisms for consideration in future policy decisions, highlighting the unpredictable nature of emerging variants and potential shifts in infection demographics associated with school closures.

Fig 1. Relative impact of school closures on SARS-CoV-2 infections, deaths and peak hospital occupancy pressure.

Results are presented as relative percentage reductions over the period 1 January 2020 to 31 December 2022 in SARS-CoV-2 infections (panels A and B), COVID-19-related deaths (C, D), and peak hospital occupancy pressure (E, F). The counterfactual “schools open” scenario was used as the reference for relative reduction calculations. In panels A, C, and E, estimates are presented as medians (horizontal lines), interquartile ranges (boxes), and 95% central credible intervals (vertical lines), and countries are listed in descending order from left to right, based on the estimated median effect for each disease indicator. The maps in panels B, D, and F present median estimates of relative percentage reductions, with negative values indicating configurations where school closures are estimated to have had an adverse impact on the considered indicator (based on the median estimate). Panel G presents the list of analysed countries along with their corresponding ISO3 codes. The maps were generated with plotly (v5.14.1) using embedded geometric data derived from the Natural Earth data set.

Fig 2. Relationship between the estimated effect of school closures on COVID-19 deaths and the total duration of school closures.

In Panel A, relative reductions of COVID-19 mortality over the period 1 January 2020 to 31 December 2022 are presented as medians (dots) and interquartile ranges (bars), with colours indicating the continents to which the countries belong. The total number of closure weeks was obtained by adding the number of weeks fully closed, and the number of weeks partially closed multiplied by 0.3 (our midpoint estimate for assumed attendance fraction during partial closures). Panel B presents the list of analysed countries along with their corresponding ISO3 codes.

Fig 3. Detailed outputs for Morocco, as a representative country where closures significantly reduced COVID-19 infections, hospitalisations, and deaths.

(A) Modelled COVID-19 deaths over time against observations. (B, C) Scenario comparison over time for COVID-19 deaths (B) and proportion infected (C). (D) Estimated transmission adjustment over time (random process) compared with school closure timelines. (E, F) Scenario comparison for cumulative deaths (E) and peak COVID-19–related hospital occupancy pressure (F) over the period 1 January 2020 to 31 December 2022. (G, H) Cumulative SARS-CoV-2 infections by age (G, age expressed in years) and strain (H) over the same period. Uncertainty is represented with dark shaded areas (IQR, interquartile ranges) in Panels A–D, and light shaded areas (95% central credible intervals) in Panel A. Boxplots E and F present estimates as medians (horizontal lines), interquartile ranges (boxes), and 95% central credible intervals (vertical lines). School closure statuses are illustrated with the horizontal coloured bands in Panels B–D (dark grey: fully closed for COVID-19, light grey: partially closed due to COVID-19, purple: academic break). In Panels A–C, the vertical lines indicate the emergence dates of the Delta variant (dashed line) and the Omicron variant (dotted line) in the country according to the GISAID database. K, thousands; M, millions; Jan, January.

Fig 4. Detailed outputs for Indonesia, as a representative country where closures may have increased COVID-19 mortality due to an exacerbated Delta wave.

(A) Modelled COVID-19 deaths over time against observations. (B, C) Scenario comparison over time for COVID-19 deaths (B) and proportion infected (C). (D) Estimated transmission adjustment over time (random process) compared with school closure timelines. (E, F) Scenario comparison for cumulative deaths (E) and peak COVID-19–related hospital occupancy pressure (F) over the period 1 January 2020 to 31 December 2022. (G, H) Cumulative SARS-CoV-2 infections by age (G, age expressed in years) and strain (H) over the same period. Uncertainty is represented with dark shaded areas (IQR, interquartile ranges) in Panels A–D and light shaded areas (95% central credible intervals) in Panel A. Boxplots E and F present estimates as medians (horizontal lines), interquartile ranges (boxes), and 95% central credible intervals (vertical lines). School closure statuses are illustrated with the horizontal coloured bands in Panels B–D (dark grey: fully closed for COVID-19, light grey: partially closed due to COVID-19, purple: academic break). In Panels A–C, the vertical lines indicate the emergence dates of the Delta variant (dashed line) and the Omicron variant (dotted line) in the country according to the GISAID database. K, thousands; M, millions; Jan, January.

Fig 5. Detailed outputs for the United Kingdom, as a representative country where closures may have increased COVID-19 mortality due to a shift in the age distribution of infections.

(A) Modelled COVID-19 deaths over time against observations. (B, C) Scenario comparison over time for COVID-19 deaths (B) and proportion infected (C). (D) Estimated transmission adjustment over time (random process) compared with school closure timelines. (E, F) Scenario comparison for cumulative deaths (E) and peak COVID-19–related hospital occupancy pressure (F) over the period 1 January 2020 to 31 December 2022. (G, H) Cumulative SARS-CoV-2 infections by age (G, age expressed in years) and strain (H) over the same period. Uncertainty is represented with dark shaded areas (IQR, interquartile ranges) in Panels A–D, and light shaded areas (95% central credible intervals) in Panel A. Boxplots E and F present estimates as medians (horizontal lines), interquartile ranges (boxes), and 95% central credible intervals (vertical lines). School closure statuses are illustrated with the horizontal coloured bands in Panels B–D (dark grey: fully closed for COVID-19, light grey: partially closed due to COVID-19, purple: academic break). In Panels A–C, the vertical lines indicate the emergence dates of the Delta variant (dashed line) and the Omicron variant (dotted line) in the country according to the GISAID database. K, thousands; M, millions; Jan, January.