Improving vaccination coverage and offering vaccine to all school-age children allowed uninterrupted in-person schooling in King County, WA: Modeling analysis

The rapid spread of highly transmissible SARS-CoV-2 variants combined with slowing pace of vaccination in Fall 2021 created uncertainty around the future trajectory of the epidemic in King County, Washington, USA. We analyzed the benefits of offering vaccination to children ages 5–11 and expanding the overall vaccination coverage using mathematical modeling. We adapted a mathematical model of SARS-CoV-2 transmission, calibrated to data from King County, Washington, to simulate scenarios of vaccinating children aged 5–11 with different starting dates and different proportions of physical interactions (PPI) in schools being restored. Dynamic social distancing was implemented in response to changes in weekly hospitalizations. Reduction of hospitalizations and estimated time under additional social distancing measures are reported over the 2021–2022 school year. In the scenario with 85% vaccination coverage of 12+ year-olds, offering early vaccination to children aged 5–11 with 75% PPI was predicted to prevent 756 (median, IQR 301–1434) hospitalizations cutting youth hospitalizations in half compared to no vaccination and largely reducing the need for additional social distancing measures over the school year. If, in addition, 90% overall vaccination coverage was reached, 60% of remaining hospitalizations would be averted and the need for increased social distancing would almost certainly be avoided. Our work suggests that uninterrupted in-person schooling in King County was partly possible because reasonable precaution measures were taken at schools to reduce infectious contacts. Rapid vaccination of all school-aged children provides meaningful reduction of the COVID-19 health burden over this school year but only if implemented early. It remains critical to vaccinate as many people as possible to limit the morbidity and mortality associated with future epidemic waves.

with states listed in Table S1 and parameters in Table S2, (1) A schematic cartoon of the model is provided as Fig S1. In words, this model describes a population of susceptible individuals S(a, v = v 0 , q), w ho can become vaccinated S(a, v ∈ v 1 , . .., v p , q) and/or infected. Infected individuals begin in an exposed state E, and after a w aiting time γ 1 proceed to infection as asymptomatic (A1, probability 1 − π) or symptomatic (P, probability π). Note that P is a pre-symptomatic period with waiting time γ 2 for symptomatic infections. After the presymptomatic period, symptomatic infections are divided into mild (IM, probability m(a)) or severe (IS, probability 1 − m(a)) infections. Mild infections recover with rate r M . Asymptomatic infections have an analogous state A2 from which they recover at rate r A . Severe infections are defined as those that lead to hospitalization, H, and possibly death, F. Hospitalized infected recover at rate r H , or die with rate f(a) and it is also possible to die outside the hospital with the same rate.
Infected individuals are diagnosed as COVID-19 cases with a rate ∆ X (a) that varies by symptom severity (∆ IM (a), ∆ IS (a), ∆ H for mild, severe, and hospitalized respectively. Asymptomatic infections are assumed to be diagnosed at a rate proportionally inferior to symptomatic defined by ρ A . Diagnosed individuals move to a set of parallel states, DA1, DA2, DP, DM, DS, DH and DF. We assume any individual who dies was diagnosed ultimately. That is, the states F and DF are combined for the total number of deaths. Diagnosed cases are used for comparing to data and for lowering transmission among diagnosed (see below).
Recovery is not a separate state, but incorporated into the susceptible state to allow for reinfection (immunity can be parameterized to be partial or complete). The number of recovered states is also flexible. F or example, information on infecting strain can be retained t o investigate cross-strain immunity, or vaccination status can be assumed to trump infection in terms of re-infection, or other possibilities. How individuals move from the infected states to the compartments of S are described by the Y matrix. Waning immunity, that is moving individuals from the vaccinated and recovered compartments of S to be newly susceptible, is also possible and is described by the W matrix.
Vaccinations occur by adjusting the status of susceptible (or recovered, because recovered individuals can still be vaccinated) within S. Infected individuals with new variants can enter the system with rate ν import (a, q). Note we assume a closed system (conservation of number), meaning this model could be interpreted as some KC resident visited elsewhere, contracted the variant, and returned or also that some KC resident left permanently precisely when a newly infected individual entered. New strains are also defined by their infectivity ν inf ( q) and severity ν sev (a, q), relative to the original strain. It is also possible to set the starting prevalence of different strains in the initial conditions, which we favor as this is generally more straightforward to estimate than the importation rate.
There are several additional model features that adjust the dynamical system, which are detailed in the following sections.

Force of infection.
The force of infection λ(a, v, q) governing transmission depends on the state of the transmitting individual (e.g. asymptomatic transmission is less likely, certain variants are more infectious, certain ages interact less). It also depends on vaccination and the time-dependent reduction in contacts mediated by social/physical distancing σ(a).
where X = {A1, A2, P, IM, IS, H, DA1, DA2, DP, DM, DS, DH} is the set of all potentially infectious states. Naturally, susceptible, exposed, recovered, and deceased individuals also do not contribute to ongoing infection.The key parameter for infection is β * , the base rate of transmission, which is further modified by state-specific transmission rates, β A 1, β A 2, β P , β M , β S , β H . Note it is assumed that hospitalized individuals do not contribute to transmission (β H = 0). For the diagnosed states, the effect of reducing interactions upon being diagnosed is handled by multiplying the statespecific β X by the reduction in transmission due to diagnosis, β D . Both β * and β D are calibrated parameters.
The model uses an empirically derived contact matrix C(a) that parameterizes the probability of interactions between transmitting and exposed individuals in different age groups. The contact matrix is assembled from those specifying interactions from different locations (Supplementary Table 3), and a location-specific social distancing parameter can also be specified. Note that these locationspecific social distancing parameters are scalars and not age-stratified.
Transmission is further affected by κ(a) the relative age-specific susceptibility. The reduction in transmission due to social distancing is handled by the σ(a), which varies between 0 and 1 and is also age-specific (i.e., the more susceptible oldest age group can have a higher value for social distancing than younger age groups). It is assumed that social distancing affects both transmission and susceptibility, thus it is applied to both sides of the contact matrix C(a). The location-specific social distancing parameters (σ school , σ home , σ work , and σ other ) reduce only interactions corresponding to that location in the contact matrix, for example due to the closing of schools.
Dynamic social distancing. We include a time-varying, age-stratified vector σ(a, t) that governs social distancing (non-pharmaceutical interventions) including reduced contacts through personal choices and/or mandated partial lockdowns, as well as reductions in exposure contacts due to mask wearing, physical distancing, capacity limitations, vaccination requirements, etc. σ(a, t) varies from 0, indicating pre-pandemic levels of societal interactivity and no masking, to 1, indicating complete lockdown with no interactions. σ(a, t), along with reductions in the susceptible population due to infection or vaccination, is the main driver of R ef f ective and therefore epidemic peaks and declines.
After the calibration period, the issue is how to set σ(a, t) for forward simulation into the future, as σ(a, t) is changing throughtime in response to humanbehavior. We use dynamical social distancing based on a specified threshold, such as the current diagnosed cases or hospitalizations or some combination. The threshold is a flexible, user-specified function, and thus can be based on government criteria for implementing restrictions, such as weekly or bi-weekly case or hospitialization numbers. This function is parameterized by the following values: the maximum T max and minimum T min thresholds for increasing or decreasing restrictions (corresponding to changing σ), the period τ over which the threshold is calculated, the restricted and released social distancing values σ max (a) and σ min , and the incrementσ inc by which to change σ when releasing restrictions. Thus we have where the sum of the threshold metric(s) τ T is taken over time period τ . Based on local policy for decision making, this is set at τ = 1 week intervals in the current simulation. Thus, for example, the system triggers increased restrictions if the weekly number of hospitalizations rises over the max threshold and distancing immediately becomes σ max (a) which is age-variable: 70% of prepandemic levels in non-seniors and 50% in seniors. Then, once hospitalizations drop below the release threshold T min , 10% of the distancing is removed every τ weeks until reaching the minimum social distancing σ min . This value is not necessarily zero because we expect persistent features such as masking, work from home and avoidance of large social gatherings will continue to limit the number of interpersonal contacts relative to pre-pandemiclevels.
Vaccination mechanisms. Original COVID-19 vaccine efficacy trials measured reductions in symptomatic disease. Therefore, it is unclear whether reductions in disease were mediated by totally prevented infections, or rather infections that were more likely to be asymptomatic. Therefore, we include the possibility of vaccines that work by several mechanisms and include three parameters controlling the effect of vaccination, VE SUSC , VE SYMP , VE INF . They are all vectors, so that if there are multiple vaccines defined, they can all have unique vaccine effectivenesses, likewise for any recovered classes. Vaccine effectivenesses are also strain-specific. The vaccine can completely block infection and reduce the number of vaccinated individuals that are susceptible by some fraction (VE SUSC ), and thus modifies the left-hand side of the contact matrix, affecting susceptibility. Or, it can block symptomatic disease in individuals who are infected despite vaccination (VE SYMP ), altering p which controls the proportion of symptomatic and asymptomatic infections. Finally, it can decrease the possibility of onward transmission in individuals who are infected despite vaccination (VE INF ), and thus it modifies the right-handside of the contact matrix, altering β * , thereby reducing transmission. Each vaccine efficacy ranges from 0-1. Furthermore, there is an additional vaccine eficacy, that against hospitalization (i.e. severe disease), conditional on being infected, VE H .
Vaccination rollout. The key parameters that govern vaccination distribution are the vaccination rate, V rate , the vaccination distribution, V dist , the vaccination coverage limit, V coverage , and the vaccination priority, V priority . The vaccination rate is the numberof vaccines distributed per day and is set per vaccine. The vaccination distribution describes what percentage of vaccines are distributed to each age group. This can be set proportional to the percentage of the population for an equal distribution or can be used to prioritize certain age groups. For example, the initial vaccine rollout was targeted primarily at the elderly, then adults, and finally certain segments of the youngest age group. The vaccination coverage describes the upper limit of coverage to apply. It is applied per age group. The vaccination priority parameter is a vector of age groups, and describes the order to reallocate available vaccines if one or more of the age groups has reached the coverage limit. All of these parameters can change throughout the vaccine rollout, i.e. to increase or decrease the doses available, using the standard temporal parameter framework.
Both susceptible and recovered individuals are eligible for vaccination, and are applied proportion-ally. While in theory individuals in the infected classes could be vaccinated, this is not implemented in the model, and this follows advice to not be vaccinated while currently ill with COVID. The V matrix describes which states are eligible for vaccination, as there can be a variable number of recovered states defined for the model, as well as what proportion of the total V rate doses available are for each vaccine.

King County model calibration.
Data. The model outputs were calibrated to three corresponding cumaltive metrics in King County WA: diagnosed cases, hospital admissions, and deaths. Thus corresponding to the cummlative values for the following states cases = (DA1 + DA2 + DP + DM + DS + DH), hosptializations = (H + DH), deaths = (F + DF). Each metric was tracked by age, which we consolidated into the 4 age groups a, so that calibration was performed against 12 metrics total.
Due to effects of weekends and weekdays, some noise in the data, and the tendency of the daily time series to be auto-correlated, we took a 7 day smoothed average of 3 days before and after the day of interest, and then used weekly values of the metrics to calibrate against (that is 1 value every 7 days). Because the metrics were on very different scales, for example with many more cases than deaths or hospitalizations and differences across age classes, we normalized the metrics. For each of the 12 metric time series, we divided by its mean, and applied the same procedure to the model output, also normalizing by the mean of the data metrics.
Algorithm. We used an Approximate Bayesian Computing (ABC) algorithm implemented in the EasyABC R package using the ABC rejection function. This algorithm samples from the prior parameter distributions, runs the model to calculate the outputs, then calculates the distance between the model outputs and the data metrics. Because we already normalized the metrics, we used a simple Euclidean distance to select the 100 best fitting parameter sets. These formed the parameter posterior distributions, and we used this set of parameterizations in the simulations to incorporate parameter uncertainty.
Calibrated parameters. The m odel calibration w as performed in two stages. An initial calibration w as performed for the initial epidemic period, w here the d ate o f the initial case t 0 w as a calibrated parameter, and the m odel w as calibrated to data through April 30, 2020. The following parameters w ere calibrated: β * , β D , t 0 , σ(a), ∆ IM ( a), ∆ IS ( a), ρ A (a), m(a). A ll parameters w ere given uniform prior distributions based on the best available data (either King County W A data or published studies). W e obtained 100 parameterizations which we used for the posterior distributions of the parameters. We performed a second calibration for the epidemic state just before the simulation period, using data from April 1, 2021 through June 15, 2021. In order to retain the initial period parameterizations for β * , β D , and m(a), we used the emprical posterior distribtion from the inital period calibration as the prior distribution for the final period c alibration. We performed the final period calibration with uniform priors for the following parameters: σ(a)[1], σ(a) [2], ∆ IM (a), ∆ IS (a), ρ A (a), where σ(a)[1] was for the period April 1, 2021 through April 24, 2021, and σ(a) [2] was for the period April 25, 2021 through June 15, 2021, in order to fit the i ncrease a nd decline in cases during the calibration period. During both calibrations the σ school parameter was set to 1, reflecting the fact that schools were closed. Results presented in Fig S2, S3, S4 and S5.

Initial conditions.
Cumulative incidence by age (ci a ) and current prevalence on June 1st, 2021 are back-calculated (see next section) from hospitalization in King County through July 31st, 2021 vaccination frequency by age group (fv a ) is derived from the CDC. Using this information, we populate all compartments using the following steps.
1. We subdivide the entire population, by age group (N a ) into three immune classes: "susceptible" (S a ), "vaccinated" (V a ), and"recovered" (R a ). This gives a total of 12 subpopulations.
2. The total prevalence by variant is set to the back calculated value for overall prevalence. It is assumed that, of initial infections, 80% are with the Alpha variant and 20% are with Delta. For each we create a 12x12 next generation matrix and use the normalized eigenvectors to subdivide the total prevalence by age and immune classes. This gives a total of 24 infected populations.
3. Each of the 24 infected populations is further subdivided into disease stages (such as symptomatic, diagnosed, hospitalized, etc.) using a quasi-steady state approximation.
All other initial populations, such as those corresponding to "waned" immunestates, are assumed to be zero initially.
Back-calculating infections . Calculation of our initial conditionsrelies on back calculation of infections from hospitalization data, which we assume to be more consistently reported. We reconstruct infections by age, x a , from the combination of hospitalizations and non-hospitalized deaths by age, y a . We assume that the time from infection to hospitalization is exponentially distributed and furthermore we penalize large differences in x a from one day to the next. We calculate y a by solving the following regularized system with non-negative least squares.

A) B)
C)