Abstract
Promising COVID-19 vaccines are under rapid development, with initial deployment expected by 2021. Careful design of a vaccine prioritization strategy across socio-demographic groups is a crucial public policy challenge given that (1) vaccine supply will be highly constrained for the first several months of the vaccination campaign (2) there are stark differences in transmission and severity of impacts from SARS-CoV-2 across groups, and (3) SARS-CoV-2 differs markedly from previous pandemic diseases like influenza. We assess the optimal allocation of a limited vaccine supply in the U.S. across groups differentiated by age and also essential worker status, which constrains opportunities for social distancing. We model transmission dynamics using a compartmental model parameterized to capture current understanding of the epidemiological characteristics of COVID-19, including key sources of group heterogeneity (susceptibility, severity, and contact rates). We investigate trade-offs between three alternative policy objectives: minimizing infections, years of life lost, or deaths. Moreover, we model a dynamic strategy that responds to changes in the population epidemiological status. Because contacts are concentrated within age groups, there is diminishing marginal returns as vaccination coverage increases in a given group, increasing the group’s protective immunity against infection and mortality. We find that optimal prioritization typically targets older essential workers first. However, depending on the policy objective, younger essential workers are prioritized to control spread or seniors to control mortality. When the objective is minimizing deaths, relative to a non-targeted approach, vaccine prioritization averts 15,000 deaths in our baseline model, with a range of 7,000-37,000 across alternative scenarios.
1 Introduction
As the novel coronavirus (SARS-CoV-2) continues to spread in many countries despite intervention efforts, public health experts see a vaccine as essential to dramatically reduce the mortality burden and possibly halt local transmission (1). Novel coronavirus disease 2019 (COVID-19) has resulted in over 930,000 confirmed deaths globally (2) as of mid-September 2020. Multiple promising vaccines are under rapid development, with deployment possible in late 2020 or early 2021 (3). While the vaccine supply schedule remains uncertain, it is clear that vaccine availability will be highly constrained for at least several months after launching the vaccination campaign (4). This scarcity, combined with stark differences in the spread and impact of SARS-CoV-2 across demographic groups, means that prioritization of the vaccine is an imminent and crucial public health challenge, and as such under active discussion by the “Advisory Committee on Immunization Practices (ACIP) of the US Centers for Disease Control and Prevention (CDC) and the National Academy of Medicine (NAM), as well as globally at the World Health Organization (WHO) and elsewhere” (5).
An effective public health policy for pandemic vaccine allocation requires an understanding of transmission and epidemiological characteristics of the novel disease across different socio-demographic groups and estimates of prevalence and built-up immunity levels when immunization begins. These key components are then integrated into a mathematical and statistical modeling framework of the state and the transmission dynamics of the novel pathogen. This framework can then be utilized to investigate the optimal vaccine allocation strategies to achieve a defined public health objective while taking into account alternative possible scenarios for uncertain components, e.g. the level of vaccine efficacy or temporal changes in vaccine availability (6, 7).
Previous experience with vaccine development mid-pandemic offers limited insights for SARS-CoV-2 prioritization. SARS and Zika vaccine development was incomplete when those outbreaks ended (8). In 2009, as the novel A/H1N1 influenza virus continued to spread across the U.S., researchers investigated optimal vaccination strategies using an age-structured dynamical model. They found that school-aged children and their parents should be prioritized, a strategy that would indirectly protect individuals at higher risk of severe health outcomes (6). Sharp differences in the epidemiology of human influenza and COVID-19 indicate that vaccination strategies against the ongoing pandemic should not simply mirror vaccination policies against influenza. For example, COVID-19 is associated with lower susceptibility to infection among children and adolescents (9) and has a substantially higher infection fatality rate overall that also increases markedly with age (10). Toner et al. (5, p. 24) provide a detailed overview of the 2018 pandemic influenza vaccination plan and conclude that, “the priority scheme envisioned…does not comport with the realities of the COVID-19 pandemic and new guidance is needed.”
We develop and apply an analytic framework to assess the optimal and gradual allocation of limited COVID-19 vaccine supply in the U.S. across socio-demographic groups differentiated by age and essential worker status. Here we use the concise label of “essential workers” to indicate workers who are (A) involved in activities essential to the maintenance of critical services and infrastructure and (B) required to work in person. These individuals comprise almost half of the U.S. workforce (see Methods). Toner et al. (5) emphasize the importance of considering the prioritization of essential workers, a group that has “been overlooked in previous allocation schemes”.
The transmission dynamics are modeled using a compartmental model tracking nine disease states as shown in Fig. 1. The parameters are set to capture our current understanding of the epidemiology of COVID-19 (see Methods). We investigate three alternative policy objectives: minimizing expected symptomatic infections, years of life lost, or deaths. We use stochastic nonlinear programming techniques to solve for dynamic vaccine prioritization policies that respond to changes in the epidemiological status of the population (shares of the population in different disease states), by updating the prioritization every month for six months. A central constraint is the scarcity of the ongoing supply of vaccines (e.g., sufficient to vaccinate 60% of the population in the first six months). We focus on the challenge of allocating vaccines across the general population, specifically the distribution of vaccine remaining after specialized top priority groups are covered (e.g. medical personnel treating COVID-19 patients).
Movement of individuals between disease states: suceptible (S), protected by the vaccine (P), vaccinated but not protected (F), exposed (E) presymptomatic (Ipresym), asymptomatic (Iasym), symptomatic (Isym), recovered (R) and deceased (D)
We account for the interaction between vaccine allocation policy and non-pharmaceutical interventions (NPI) by modeling location specific transmission rates reflecting social distancing patterns caused by increases in remote work, school closures and reductions in social engagements (see Methods). A key component is the scaling of transmission in the workplace to reflect the different constraints faced by essential versus non-essential workers in their capacity to engage in social distancing. We also scale the probability of transmission from any given contact to capture the effect of non-social-distancing (NSD) NPI such as the use of masks, hand washing and maintaining physical distance while meeting. We construct a Base scenario to establish a set of plausible outcomes and then analyze a range of alternatives to explore key uncertainties about the state of the world upon a vaccine’s arrival. We consider variation in the following: the strength of NSD NPI; the overall and age-specific efficacy of the vaccine; the susceptibility of younger individuals to infection; the available supply of vaccine; the social contact rate among school-age children; and the mechanism by which vaccines work (by reducing transmission and/or symptoms).
To our knowledge there are no published analyses of optimal COVID-19 vaccination prioritization. Analysis in preprint form includes Matrajt et al. (11) and Bubar et al. (12).1 Both consider the optimal allocation of vaccines across five or more age groups. Their approaches feature rich exploration of policy sensitivity to vaccine efficacy and availability. Matrajt et al. is particularly detailed in this respect, while Bubar et al. uniquely consider differences in demographics and contact rates across multiple countries. Our analysis is differentiated by a deeper approach to the behavioral, demographic and decision models by addressing social distancing, essential worker groups, and allocation policies that can change over time.
General ethical guiding frameworks for vaccine prioritization decision-making have appeared earlier in the literature. Toner et al. (5) emphasize promoting three ethical values: the common good; fairness and equity; and legitimacy, trust and communal contributions to decision-making. Emanuel et al. (4) promote four ethical values: maximizing benefits, treating equally, instrumental value, and priority to the worst off. Our analytic focus on minimizing new infections, years of life lost, or deaths emerges from promoting “the common good” or “maximizing benefits”. Our focus on essential worker groups illustrates how ethical values (e.g. prioritizing essential workers due to the fairness of protecting those placing themselves at risk) may overlap with the common good (e.g. prioritizing essential workers to best reduce mortality and transmission). Issues of fairness and equity and protecting the worst off are not directly analyzed here but remain critical considerations.
To capture, at least coarsely, the likely feedback between falling infections and relaxation of social distancing, we assume that this distancing partially attenuates at a given threshold for new infections (see Methods). However, for the sake of simplicity, we do not address in detail the potential set of complex and differential feedback processes between health status and opening of schools, workplaces and other institutions. While we limit policy objectives to a concise metric of health outcomes (minimizing expected cases, years of life lost, or deaths) we acknowledge that other values of returning to school, work and social life are important. Finally, we do not address additional vaccine complications, such as temporary efficacy, potential side effects or any failure to take a second dose of the vaccine if necessary.
Previewing our results, we find that optimal allocation strategies are responsive to both the initial and evolving epidemiological landscape of the disease. When deaths are considered, vaccines are initially allocated to older essential workers 40 − 59 yrs. followed by seniors 60+ yrs.; when years of life lost are minimized both ages groups of essential workers are prioritized followed by seniors 60 − 74 yrs.; and when infections are minimized essential workers are prioritized followed by school-age children. In general, we find that these results are robust across a range of possible scenarios and parameters sets. However, they are sensitive to changes in vaccine efficacy and susceptibility between age groups, indicating that these uncertainties are priorities for research emphasis prior to vaccine distribution. With respect to key outcomes, we found that the optimal strategies outperformed a non-targeted strategy (e.g., distributed proportional to the size of each group) by 8 to 20% for a given target policy objective. For example, when the focus is minimizing mortality, we find that optimal vaccine prioritization averts 15,000 deaths relative to a non-targeted approach in our baseline scenario, with a range of 7,000 to 37,000 lives depending on the effectiveness of the vaccine and concurrent non-pharmaceutical interventions.
2 Results
We present results for the Base model with essential worker demographic groups and then show the sensitivity of these results with respect to the alternative scenarios. In Fig. 2 we show the optimal allocation of vaccines given each objective for the Base model. The allocations are shown on a monthly basis for six decision periods and then cumulatively (in percent of vaccine and percent of group vaccinated). Broadly, we find that the optimal policy is very dynamic: specific groups are targeted each period and these targets shift over time. Furthermore, targeting is very narrow and becomes less so as a large fraction of the population has been covered.
The optimal allocation of vaccines (vertical axes) between demographic groups for each decision period (horizontal axis in panel A-C) under the Base scenario. The three rows represent each objective, to minimize deaths (A), minimize years of life lost (YLL, B) and minimize infections (C). The bars for the six decision periods show the percentage of vaccines allocated to a specific group (indicated by a letter, color, and stars indicating essential worker groups) in that period. The two final columns show cumulative measures at the end of month six, the percent of vaccines allocated to each group (D-F), and the percent of each group that has been vaccinated (G-I). The whiskers on each bar represent the sensitivity of the optimal solution to small deviations in the outcome, specifically the range of allocations resulting in outcomes within 0.25% of the optimal solution.
In general we find that optimal dynamic allocation does not cover 100% of the susceptible population in any single demographic group before switching to another group. Further, an allocation may initially prioritize one group, only vaccinate a fraction of that group and then prioritize that group again two or more decision periods later. These switching dynamics emerge because there are diminishing marginal returns to vaccinating additional individuals in a group as more of that group becomes vaccinated.
The whiskers on bars in Fig. 2A-C show the range of alternative allocations that still produce an outcome that is within 0.25% of the optimum. For example, in the first period of Fig. 2A the whiskers show that some limited substitution in the allocation between groups d∗, f ∗ and h (but not others) can occur without a substantial reduction in the optimized outcome (minimizing deaths). In general, we find that these whiskers become more pronounced as periods progress. This shows that it becomes less critical to precisely follow the exact optimal allocation as vaccine coverage of the population expands. We also find that despite targeting of certain groups in certain periods, after six periods the percent vaccine allocated to each group and the percent of each group vaccinated (Fig. 2D-I) is more even across most groups. Overall, we find that pre-school-age children are substantially less targeted than most groups (conditional on having relatively few contacts and lower susceptibility in the Base model).
Across objectives there are substantial differences in which groups are targeted early on. When minimizing deaths, targeting progresses from older essential workers (40-59*), to the oldest (75+), to younger seniors (60-74), and then younger essential workers (20-39*). These groups are a mix of those at high risk of mortality (older groups) and high risk of contraction and spread (essential workers). When minimizing YLL, younger essential workers and younger seniors are targeted earlier (given their longer average years of life remaining). Finally, when minimizing infections we find that younger essential workers take top priority, followed by older essential workers and school-age children (5-19), since these groups have higher contacts and thus risk of contraction and spread.
Results for the “Age-only” scenario—which does not distinguish essential workers—are broadly similar (see SI Appendix). However, a significant difference is that the essential worker formulation presented in the main text targets essential workers before other working age adults and prioritizes these groups before the higher risk, 60+ age groups.
In Fig. 3A we show the dynamic path of infections, starting from the period in which vaccines become available, under various policies. As expected, infections are highest given no vaccines. Results for allocating vaccines in a manner “proportional” to each group’s size shows the substantial value of vaccines even with no targeting. As expected, the policy for minimizing infections leads to the lowest level of infections. In Fig. 3B we show the performance of each targeting policy relative to outcomes achieved with a proportional allocation, or each objective after 240 days (the policy for the sixth month is used from days 150 to 240). Overall we find that when focusing on minimizing a particular outcome, that outcome is reduced by 17-18% (relative to a proportional allocation). In the first cluster of bars, as expected the policy that minimizes deaths (“D” in green) leads to the greatest reduction in deaths (18% fewer than under a proportional allocation). However, trade-offs are stark in certain cases: the third cluster of bars shows that minimizing deaths involves a strong opportunity cost in terms of infections, which are higher than even under the proportional policy. The YLL policy is most consistent, performing second best (and never third best) when considering other outcomes (infections and deaths).
The number of infections per 1,000 individuals over time under each policy (a) and the performance of each policy relative to a uniform allocation policy (b) for the Base scenario in the essential worker model. The bars are clustered by the resulting outcome and colored by the objective driving each policy.
2.1 Sensitivity of prioritization across scenarios
To assess how robust the findings are to key uncertainties in the model we solved for optimal allocation across a range of alternative scenarios. The differences between these scenarios and the Base case are detailed in Table 1. Iterations of Fig. 2 and Fig. 3 for all alternative scenarios are shown in SI Appendix B. To compare and contrast cumulative vaccination results, in Fig. 4 we show for each of the alternative scenarios the percentage of each group vaccinated after three months (A-C) and six months (D-F). In general we find differences across groups that lessen (but not completely) by month six. We also find some differences across scenarios and objectives that differ by the horizon considered.
Descriptions of alternative scenarios to the Base model (see SI Appendix A.3 for specific levels).
The cumulative percent of each demographic group (x-axis) vaccinated after three months (A-C) or six months (D-F) under alternative scenarios (y-axis) and each objective (A, D minimize deaths; B, E minimize YLL; and C, F minimize infections). Vaccine supply is sufficient to vaccinate 10% of the population each month, except in the “Ramp up” scenario which involves 5% per month for three months, then 15% per month.
Certain scenarios are distinctive. For example, when the lower susceptibility enjoyed by those under 20 in the Base scenario is replaced by “Even susceptibility” for all, over the first three months we see substantial substitution to school-age children and away from older essential workers (min. deaths or infections) or younger seniors (min. YLL). In a second example, both weaker vaccine scenarios do not substantially change results at three months, except if minimizing deaths, in which case vaccination shifts from younger seniors to younger essential workers. When fewer vaccines are available in the first three months due to a “Ramp up” in supply, the deficit mainly accrues to older essential workers (min. deaths or infections), younger seniors (min. deaths or YLL), younger essential workers (min. YLL), or school-age children (min. infections).
At three months, the only two groups consistently not targeted are pre-school-age children and older non-essential workers. By six months, there are two groups consistently targeted: older and younger essential workers. Across each objective, one of the essential worker groups has either the highest or second highest coverage rate. The only exception to this rule is if school-age children are equally susceptible to infection as adults (as discussed above).
2.2 Vaccines partially effective at the individual level
In our results discussed above we have assumed that vaccines are completely effective at preventing infection for a given percentage of each group (e.g. 65%). An alternative approach is to consider vaccines as partially effective for everyone at the individual level. Furthermore, it may be the case that a vaccine is more effective at reducing symptoms than preventing infections. To analyze this case we extended our model structure in Fig. 1 to separate infected individuals into those previously vaccinated or not. For vaccinated individuals, we replaced the single vaccine efficacy parameter with separate parameters for reduction in spread (scaling susceptibility and transmissibility) and mortality (scaling infection fatality rate). We consider three cases. First, for comparison with the Base model, we consider the same level of efficacy (65%) for spread and mortality. Finally, we consider two variants in which the vaccine is more effective at reducing mortality but less effective at preventing spread. Specifically we model an extreme case where the vaccine reduces susceptibility and infectiousness by 10% and infection fatality rate by 90% and a moderate case where these values are 30% and 70% respectively.
Details on modeling and results for these cases are provided in SI Appendix C. In general, results are similar between the Base model (65% of individuals 100% protected by vaccination) and the three scenarios we tested with the partially effective vaccine model formulation. However, when the vaccine had an efficacy of 65% for spread and mortality, and when the vaccine had an efficacy of 90% against mortality but only 10% against spread, the optimal solutions allocated a greater share of vaccines to ages 60+ to minimize years of life lost and deaths compared to the Base model. Surprisingly, the scenario with a vaccine that reduces mortality by 70% and susceptibility and infectiousness by 30% was qualitatively different to the others, shifting prioritization towards younger essential workers and away from older non-essential workers. This finding illustrates that the indirect benefits of reducing spread (e.g. in younger essential workers) is still important to take into account even if the vaccine is relatively less effective at mitigating infections than deaths.
3 Discussion
Key insights and results from our analysis are summarized in Box 1. Together these lessons show the strong implications of considering dynamic solutions, social distancing and essential workers (given their limitations in social distancing) for vaccine prioritization.
While vaccine prioritization discussion often takes the form of identifying tiers that should be vaccinated to completion before moving on (e.g., see (5, p. 25)), we find that the optimal approach does not involve seeking 100% coverage in a single group before prioritizing other groups. In fact, in some cases a group is prioritized early on, and then revisited two or more periods later. These findings are indicative of the diminishing marginal returns to vaccinating individuals within a demographic group. Because social contacts are concentrated within groups, as vaccine coverage in a given group increases, risks for that group fall more sharply than for others. Thus, as group vaccine coverage increases, we see attenuation in the direct benefits (protection of the vaccinated) and indirect benefits (protection of the unvaccinated) of additional vaccinations in that group.
Key insights and results
Benefits: Prioritization can reduce a particular undesirable outcome (deaths, YLL, or infections), by 17-18% in the Base scenario (or 8-23% depending on the alternative scenario).
Dynamic prioritization: (A) Prioritization is responsive to the initial and evolving disease status; (B) Diminishing marginal returns to additional vaccination within a group drives a shift to other groups well before 100% vaccination of the first group is achieved.
High prioritization: Under the Base scenario, group prioritization over the first three months starts with essential workers then, depending on the objective, progresses to older and younger seniors (deaths), younger seniors only (YLL) or school-age children (infections).
Low prioritization: At three months (30% of the population vaccinated) the only two groups consistently not targeted are pre-school-age children and older non-essential workers.
Widening prioritization: As vaccination rates rise, precise prioritization becomes less critical and targeting widens to a larger set of groups.
Sensitivity: Some prioritization results are sensitive to the scenario modeled—especially whether or not those under 20 are less susceptible.
Trade-offs: Focusing on one objective leads to sacrifices in the other objectives, typically strongest when minimizing infections.
Available existing analysis of optimal COVID-19 vaccination targeting in preprint form is limited to Matrajt et al. (11) and Bubar et al. (12). Before comparing and contrasting results some key modeling differences should be noted. Both preprints consider a wider range of vaccine availability than considered here. Our analysis uniquely incorporates non-pharmaceutical interventions (NPI), including social distancing and non-social distancing (e.g. mask wearing). Doing so allows us to account for differences between groups like essential workers constrained in distancing versus others who are much less so. Matrajt et al. and Bubar et al. both implement static optimization where all vaccination available is allocated and administered in a one-shot process. Our allocation is dynamic, responding to changing conditions over a six-month period. Finally, Matrajt et al. and Bubar et al. model vaccines as “leaky”, i.e., reducing the probability that a susceptible individual will be infected. Bubar et al. also considers an “all-or-nothing” vaccine that is 100% effective for a fraction of the population. In our Base model the vaccine is “all-or-nothing”, though we also consider a leaky vaccine, as discussed at the end of the Results.
Matrajt et al. (11) found that optimal strategies to minimize deaths and years of life lost will either exclusively target groups with high infection fatality rates maximizing the direct benefit of vaccines, or will target groups with high rates of infection maximizing the indirect benefits of the vaccine. In contrast, our results indicate that optimal policies initially target groups with high risk of infection and then switch to targeting groups with high infection fatality. This difference most likely follows from our dynamic versus static allocation. The switching behavior we identify is consistent with past work on pandemic influenza vaccine prioritization, which suggests that early in an outbreak when the infection rate is growing targeting spread (maximizing indirect benefits) is more efficient, but later when the infection rate is leveling off or declining maximizing direct protection is most efficient (18). Consistent with this explanation, we find that the oldest group is prioritized in the first decision period when deaths are minimized in the “Strong NPI” scenario where the number of infections are declining, compared to the Base case where they are not prioritized until the second period.
Bubar et al. found that prioritizing adults older than 60 years of age is a robust strategy for minimizing deaths. In contrast we find that working-age adults are a key priority group, particularly older essential workers. These differences may either arise from our allowance for social distancing and/or dynamic allocation. Our accounting for social distancing on COVID-19 transmission increases the modeled benefits of targeting essential workers, who are less able to substantially reduce their social contacts than individuals over 60. Furthermore, as discussed above, the ability of dynamic policies to switch over time allows the allocation schemes we discuss to capture the benefits of using the initial vaccine supply to slow transmission without sacrificing direct protection of more vulnerable individuals later on.
While we explored a large set of alternative scenarios, there are other important possibilities that we have not included. For example, if certain population groups (e.g., children or seniors) experience significant side effects from the vaccine, prioritization might shift away from these groups (19). Another key component is the set of logistical constraints imposed by the distribution network used. Vaccines will likely be administered through various points of contact with the community (pharmacies, clinics, schools, etc.). For some demographic groups there may be differences between the share of vaccines targeted to that group and the actual share received, e.g., due to constraints in prediction and implementation.
We do not consider the potential for vaccine hesitancy in the model. In general, we find that it is not necessary or even ideal to vaccinate all of the susceptible individuals in a demographic group, at least given the limit of 60% of the population vaccinated we considered here. Thus, at least initially, some level of vaccine hesitancy may not have a material impact. However, hesitancy may play a more significant role in the longer run, especially if hesitancy rates are large and herd immunity proves difficult to achieve (e.g. if vaccine efficacy is low, and/or NPI relaxation is aggressive). Vaccine hesitancy that is concentrated in a particular community or demographic group could also change the optimal prioritization strategy. Similarly, adjustments would be needed if groups differ in the duration of vaccine efficacy or diligence in obtaining a second dose of the vaccine (if necessary).
For simplicity we limited policy objectives to a set of concise metrics of health outcomes (minimizing expected cases, years of life lost, or deaths). However, other health-related metrics such as protecting the most vulnerable and social values such as returning to school, work and social life are important to consider. Our analysis reveals that optimal strategies for minimizing deaths and years of life lost are broadly aligned with the goal of protecting the most vulnerable. These solutions target essential workers who are the least able to participate in NPI such as social distancing and thus are the most at risk of infection, and individuals over the age of 60 who have the highest risk of deaths if infected by the disease. Other social values such as returning to school will most likely change the allocation schemes to offset the risk created by relaxing social distancing. For example, if allowing children to return to school was a high priority, then allocation strategies might be tilted towards targeting school-age children and teachers. A detailed analysis of optimal vaccine allocation given the relaxation of social distancing to achieve particular social objectives is a promising direction for future research.
4 Methods
4.1 Model
To investigate the impact of vaccination strategies on the COVID-19 pandemic in the U.S., we employed a structured compartmental transmission model similar to (20). We incorporated the demographic structure of the population by tracking six age groups in the set J = {0-4, 5-19, 20-39, 40-59, 60-74, 75+}. We then extend this set to differentiate essential workers by splitting the two prime working age groups into two groups—non-essential workers (20-39, 40-59) and essential workers (20-39*, 40-59*)—yielding four groups of prime working age individuals and a total of eight demographic groups in J = {0-4, 5-19, 20-39, 20-39∗, 40-59, 40-59∗, 60-74, 75+}. For each demographic group we tracked 9 epidemiological states: susceptible (S), protected by a vaccine (P), vaccinated but unprotected (F), exposed (E), pre-symptomatic (Ipre), symptomatic (Isym), asymptomatic (Iasym), recovered (R) and deceased (D). In Fig. 1 we display the compartmental diagram and directions of transitions between epidemiological states.
We modeled the COVID-19 transmission dynamics using a system of coupled ordinary differential equations for each demographic group, indexed by i and j. The transmission rate was given by the product of the transmission probability (q), the age-specific susceptibility (si), strength of non-pharmaceutical interventions (θ), the relative infectiousness of each symptom type (τm)—where m ∈ M ≡ {asym, pre, sym}—and the rate of contact (rm,i,j) between infected individuals with symptom type m from group j and susceptible individuals from group i. The exogenously given population vaccination rate at time t is given by v(t), where units of time are days.2 In our Base model we assume that for each individual the vaccine either works or it does not (though we also consider vaccines that are partially effective for all vaccinated in our sensitivity analysis). Individuals in group i are vaccinated at a rate of µiv(t) and a fraction of the those (ϵi) are protected while a fraction remain susceptible and move to the failed vaccination category (F).3 Once infected, individuals move from exposed to pre-symptomatic at rate 
. Pre-symptomatic individuals become symptomatic or asymptomatic at rates σasym/Dpre and (1 − σasym)/Dpre, respectively. Asymptomatic individuals recover at an uniform rate 
and symptomatic individuals either recover or die at a rate of (1 − δa)/Dsym or δa/Dsym, respectively, where δa is the age-specific infection fatality rate. These assumptions yield the system of differential equations described in SI Appendix A.1, with parameter values given in SI Appendix A.2 and A.3.
4.2 Contact rates
Contact rates indicating the level of direct interaction of individuals within and between groups drive the transmission dynamics in the model. We built the contact matrices used in this model from the contact matrices estimated for the U.S. in (21). These estimates are given for age groups with five year age increments from 0 to 80 yrs. These estimates were aggregated to provide estimates for the coarser age structure used in our model. We also extended these data to estimate the contact rates of essential workers. A detailed derivation of these contact rates can be found in SI Appendix A.6 but, in short, we assumed that essential workers have on average the same pattern of contacts as an average worker in the population in the absence of social distancing. We then scaled the contact rates for essential and non-essential workers to represent the effects of social distancing and calculated the resulting mixing patterns assuming homogeneity between these groups.
We constructed contact matrices for each of four locations, x ∈ {home, work, school, other}, following (21). The total contact rate for an asymptomatic individual before the onset of the pandemic is given by the sum of these location specific matrices. However, it is clear that populations are exhibiting social distancing in response to the pandemic (22). We further expect symptomatic individuals to change their behavior in response to the illness. We account for these behavioral changes as described next.
4.3 Social distancing
Expression of symptoms and social distancing policies are likely to change individuals’ behaviors over time. To model these changes we scaled the contribution of each contact matrix for location x:
The weights αm,x depend on disease and symptom status (m) and location (x) as specified in Table 2. We scaled social contacts for symptomatic individuals following changes in behavior observed among symptomatic individuals during the 2009 A/H1N1 pandemic (23). For those without symptoms (susceptible and asymptomatic) the weights were specified to match reduced levels of social contacts as the product of social distancing policies. Home contact rates were held constant. Since completed research studies to understand changes in work contact rates are not yet available, we select a level based on preliminary survey data across eight U.S. regions collected by the Institute for Transportation Studies at the University of California Davis, which indicates that trips to work have fallen after the onset of the pandemic from an average of 4.1 to 1.9 days, or 54% (24). The work contact rates for both model formulations were set to be consistent with an overall reduction of 54%. This value was used directly in the “Age-only” scenario and divided into a weight of 100% for essential workers and a weight of 9% for non-essential workers in the Base model (with essential workers). School contact rates were set to an assumed weight of 30%, to account for a mixed effect of a small fraction of schools remaining open and possible increased social contacts between school-age children during time that would otherwise have been devoted to school. As an alternative scenario, we considered the case of more school contacts with a weight of 70% (see SI Appendix Table 4). Contacts in other locations were given an assumed weight of 25%.
Weights on contact rates for a given disease and symptom type (m) and location/activity (x) under social distancing. When essential and non-essential-worker weights are both needed the former is marked with a star.
Parameter values that differ from the Base case for alternative scenarios.
The proportion of workers deemed essential, p, was estimated with two components: the total number of workers involved in activities essential to the maintenance of critical services and infrastructure and the fraction of these workers that were required to work in person. The cyber-security and infrastructure security agency of the U.S. estimates that 70% of the workforce is involved in these essential activities (e.g. healthcare, telecommunications, information technology systems, defense, food and agriculture, transportation and logistics, energy, water, public works and public safety)(25). We used estimates of the fraction of workers that could successfully complete their duties from home produced by (26) who estimated this value at approximately 30%. These two values gave a final proportion of p = 0.7(1 − 0.3) = 0.49.
4.4 Transmission rate calibration
The model was calibrated to match the observed R0 for COVID-19 in the U.S. (see SI Appendix Table 3) by solving for probability of transmission q, assuming a naive (pre-pandemic) population. Details of this procedure are provided in SI Appendix A.4.
Base model parameter values and sources.
4.5 Initial conditions
Because the expected epidemiological conditions {Ipre(0), Iasym(0), Isym(0), S(0)} by the time the initial vaccine doses are ready for deployment are uncertain, we constructed plausible baseline values for the U.S. using estimates of COVID-19 disease burden from the start of the outbreak in February 2020 through mid-September 2020 and used projected disease burden estimates by December 1, 2020 taken from near real-time projections by the Institute of Health Metrics and Evaluation (27). Specifically, we set the initial epidemiological conditions to be consistent with cumulative and current cases by December 1, 2020. These cases were apportioned between demographic groups to reflect the attack rates of COVID-19 for each group under the given social distancing policy.
4.6 Vaccine prioritization optimization
The planner’s decision problem is to allocate the daily supply of vaccine (v(t)) across the demographic groups according to a given objective. We assume that this group allocation vector, µ, can be chosen on a monthly basis at the beginning of each of the first six decision periods. We numerically solved for vaccine allocation strategies that minimize the total burden associated with three different health metrics: deaths, years of life lost (YLL) or symptomatic infections:
where ei is the years remaining of life expectancy for group i and with a six month time horizon (T = 180 days). Preventing deaths and years of life lost are “consensus value(s) across expert reports” (4, p. 2052) while some argue that “protecting public health during the COVID-19 pandemic requires…minimizing COVID-19 infection” (5, p. 10).
We solved for the optimal allocation of available vaccines across demographic groups for each month over six months. We identified the optimal solution using a two-step algorithm. In the first step we used a genetic algorithm similar to (28) to identify an approximate solution. This approach uses random sampling of the potential solution space to broadly explore in order to avoid narrowing to a local and not global minimum. In the second step we used simulated annealing to identify the solution with precision. At a given optimal solution, it may or may not be the case that the outcome of interest (e.g. minimizing deaths) is sensitive to small changes in the allocation decision. Thus, around the optimal allocation we also identified nearby allocations that produce outcomes that are less desirable but still within 0.25% of the optimized outcome. A detailed description of the algorithm is given in SI Appendix A.7. All code for the optimization was written in the Julia programming language (29).
Data Availability
All data used in the article are publicly available from the relevant citation included in the text. Code for the analysis may be obtained from the authors upon request.
Data access
All data used for informing the numerical analysis are freely available at the source noted for each measure. The data and code used to initialize and run the models will be made available on Github upon publication, https://github.com/XXX.
A Model specification, parameterization and optimization
A.1 Model dynamic equations
The dynamic equations specifying transitions between the disease states are as follows:
To reduce clutter we have suppressed the time index t on each of the state variables, the vaccine allocation vector µi, and the vaccination rate v.
A.2 Model parameters
A.3 Parameters for alternative scenarios
A.4 Calibration
The relationship between the basic reproduction number, R0, and parameters governing transmission and epidemio-logical characteristics is given by the so-called next-generation matrix:
where the maximum eigenvalue operator wraps several terms including r, the social contact matrix, s, the age-specific susceptibility rate, n, a vector of the proportions of the population in each demographic group and Δ, an operator that signifies multiplying each row of a matrix by the corresponding entry in the vector. For symptom type m ∈ {asym, pre, sym}, the constants Dm, τm and σm represent the duration, relative infectiousness of an individual and the probability of type m, respectively.
We first set a baseline R0 = 2.5 as estimated by (37). We then solve for the transmission probability parameter, q, using Equation 14, assuming a naive (pre-pandemic) population. We then scaled q by a fixed factor θ ∈ [0, 1] to reflect the impact of non-pharmaceutical interventions (NPI) like masks, hand washing and maintaining distance when contacts are made.
A.5 Initial conditions
We set plausible baseline values for initial epidemiological conditions that will be present in the U.S. population when a vaccine arrives using estimates of actual U.S. cases from the start of the outbreak (February 2020) through present and expected cases projected forward through December 1, 2020 produced by the Institute for Health Metrics and Evaluation (27). For the share of the population that is symptomatic, Isym(0), we take the estimated case count for the Dsym = 7 days (the duration of symptomatic infection) ending December 1, 2020, multiply by the share of cases that are symptomatic (1 − σasym) and then divide by the U.S. population. For the pre-sympt omatic share, Ipre(0), we assume it is consistent with Isym(0) but scaled by its shorter duration: 
. The asymptomatic share is calculated by scaling Isym(0) to account for differences in duration and the relative share of asymptomatic to symptomatic cases: 
. Finally, we estimate the share of the population that is recovered, R(0), by taking the cumulative projected case count through December 1, 2020, divided by the total population.
A.6 Contact matrices distinguishing essential workers
Estimated contact rates for the U.S. were obtained from (21) who used population-based contact diaries from the European POLYMOD survey to project to other countries, including the U.S. These included contact rates for 16 age classes in five year increments from ages 0 to 80. We collapsed these to five age groups (0-4, 5-19, 20-39, 40-59, 60-80) using population-weighted sums:
where {i, j} are the subscripts for the five year age bins, {i, j} are the subscripts for the larger age bins, ri,j,x is the average number of daily contacts a person in group i makes with a person in group j for activity/location x, and 
is the population size for age group h.
The total number of i-to-j contacts must equal the total number of j-to-i contacts: 
Because numerical issues—estimation in (21), bin discretization and rounding—can lead to small differences, we ensure this condition holds by imposing,
where the numerator is the mean of the two measures of total contacts between groups i and j and the denominator transforms the result to per-capita in i.
Setting essential worker contact rates requires additional assumptions and attention to the activity/location. We define the essential worker indicators e ∈ {n, y} for “no” and “yes”. Our grouping is such that all essential workers (e = y) are employed but non-essential-workers (e = n) are a mix of employed and not employed. Let e′ represent the indicator for a second group which can be equal or not equal to the value for e.
In the case of all activities/locations x that are not work, contact rates are given by
This follows from the assumption that contacts made by any group (i, e) with any other group (j, e′) are independent of i’s essential worker status. Thus, we only need to split contacts ri,j,x into those made with essential worker type e′ = y versus the remainder with type e′ = n, i.e. given the share 
.
Estimating contacts when x = work involves a larger number of steps. We first address contacts made by essential workers (e = y) before turning to non-essential workers (e = n). For e = y, let the share of the working age population (20 − 59) in group i that is employed be given by pi.
We assume that all of the work contacts are attributable to employed adults resulting in an employed adult contact rate of ri,j,work/pi. Then the contact rate of essential workers (e = y) in group i with age group j is
Let the fraction of working age group i that is employed in an essential worker role be given by pi,y. The average workplace contact rate for non-essential-workers in group i with group j is given by
where αwork < 1 scales for social distancing and the final term in brackets scales for the share of non-essential workers that are employed and thus have contacts at work.
Finally, we assume that the average work place contact rate for an individual of type (i, e) with individuals of type (j, e′) is given by the partial contact rate r(i,e),j,work times the proportion of total work contacts of individuals in group j that are made by individuals in sub group e′:
Finally we scale the work contacts for age groups that are not separated into essential and non-essential workers (5-19, 60-80) to match the scaling for prime working age classes.
A.7 Optimization algorithm
The optimization algorithm used in our analysis is split into two parts. First a genetic algorithm is run to identify an effective strategy near a global optimum. This solution is then refined via simulated annealing.
Genetic algorithms take inspiration from the natural process of evolution, and work by randomly sampling a populating of candidate solutions, selecting a set of survivors based on the candidates performance against the objective function, information from these survivors is then used to generate a new generation of candidates solutions, and so forth (28). The genetic algorithm executes the following steps:
Sample Nt=0 candidate solutions {xn,t=0} from a Dirichlet distribution with parameter α0.
Each candidate solution is evaluated with the objective function.
The bests Kt=0 candidates
are solved and the distributions parameter α0 is updated to α1 by maximizing the likelihood of of the .Steps 1 to 3 are repeated for a fixed number of iterations T and the best candidate solution sampled at any iteration is returned. The values Nt and Kt are tuned for each step to maximize performance.
Simulated annealing is based on thermodynamic models of cooling metals. Briefly, the algorithm is initialized by sampling a candidate solution x0, this candidate solution is updated by sampling a new candidate solution xt from a proposal distribution centered around x0. This solution is either accepted and replaces the current x0 or it is rejected and a new candidate solution is drawn using the existing value of x0. The proposed solutions xt are accepted if they perform better against the objective than the incumbent x0, if xt > x0 it is selected with probability p = exp [−(xt − x0)/Tt]. large values of Tt increase the probability that a new candidate solution will be accepted allowing the algorithm to explore the solution space and move away from local minima. Tt is reduced over time to allow the algorithm to start exploring the solution space and then eventually stabilize on a global minimum. The simulated annealing executes the following steps:
Initialize a chain with value x0. Generate a new sample from the proposal distribution xt ∼ tr(N (tr−1(x0), σI)) where the transform tr from Rn to the solution space. Initialize a counter i that tracks the number of iterations.
If xt < x0 replace x0 with xt, update i = i + 1 and repeat from step 1.
If xt > x0 sample µ ∼ unif (0, 1). If µ > exp (−(xt − x0)/T (i)) then replace x0 with xt update i and repeat from step 1. Otherwise save x0 and repeat from step 1. We used T (i) = T0/i as the temperature function.
Stop when i > max iter
These algorithms were tuned experimentally to consistently converge on a minimum solution on a test case. We used the minimum years of life lost under the Base parameter set as our test case.
To quantify the sensitivity of the solutions to deviations in the outcome of interest, we sampled solutions near the optimal solution using a Markov chain. This algorithm is initialized at the optimal solution identified by simulated annealing and the genetic algorithm and samples are drawn from a proposal distribution and accepted if they perform within the desired tolerance (0.25%) of the optimal solution.
B Results for additional scenarios
In this section we present detailed results for each alternative scenario. For simplicity we do not repeat extended figure captions, which follow those of main text Figures 2 and 3 (except for the stated scenario).
B.1 Age-only demographic groups (without essential workers)
The optimal allocation of vaccines between age groups at each decision time point. Each panel shows the solution for a given objective.
The number of infections per 1,000 individuals over time under each policy (left) and the performance of each policy relative to a a uniform allocation policy (right).
B.2 Strong NPI
The optimal allocation of vaccines between age groups at each decision time point. Each panel shows the solution for a given objective.
The number of infections per 1,000 individuals over time under each policy (left) and the performance of each policy relative to a a uniform allocation policy (right).
B.3 Weak NPI
The optimal allocation of vaccines between age groups at each decision time point. Each panel shows the solution for a given objective.
The number of infections per 1,000 individuals over time under each policy (left) and the performance of each policy relative to a uniform allocation policy (right).
B.4 Weak vaccine
The optimal allocation of vaccines between age groups at each decision time point. Each panel shows the solution for a given objective.
The number of infections per 1,000 individuals over time under each policy (left) and the performance of each policy relative to a uniform allocation policy (right).
B.5 Weak vaccine in seniors
The optimal allocation of vaccines between age groups at each decision time point. Each panel shows the solution for a given objective.
The number of infections per 1,000 individuals over time under each policy (left) and the performance of each policy relative to a uniform allocation policy (right).
B.6 Lower child susceptibility
The optimal allocation of vaccines between age groups at each decision time point. Each panel shows the solution for a given objective.
The number of infections per 1,000 individuals over time under each policy (left) and the performance of each policy relative to a uniform allocation policy (right).
B.7 Even susceptibility
The optimal allocation of vaccines between age groups at each decision time point. Each panel shows the solution for a given objective.
The number of infections per 1,000 individuals over time under each policy (left) and the performance of each policy relative to a uniform allocation policy (right).
B.8 Ramp up
The optimal allocation of vaccines between age groups at each decision time point. Each panel shows the solution for a given objective.
The number of infections per 1,000 individuals over time under each policy (left) and the performance of each policy relative to a a uniform allocation policy (right).
B.9 Open schools
The optimal allocation of vaccines between age groups at each decision time point. Each panel shows the solution for a given objective.
The number of infections per 1,000 individuals over time under each policy (left) and the performance of each policy relative to a a uniform allocation policy (right).
C Vaccine that is partially effective at the individual level
Vaccines can provide multiple forms of protections against infections. Among these protections is the ability for vaccines to prevent individuals from becoming infected (the case considered in the main text). In addition, if vaccinated individuals still become infected they may (1) develop reduced infectiousness and/or (2) exhibit reduced severity. To allow for these latter two cases we changed the model structure to track vaccinated and infected individuals. To do this we maintained the protected and uninfected category P and added four categories: vaccinated and exposed class Pexp, vaccinated and pre-symptomatic Ppresym, vaccinated and asymptomatic Pasym and vaccinated and symptomatic Psym. The three vaccine efficiencies are modeled by three age specific vectors, V Esucpt, V Etrans, and V Esym, which represent the amount the vaccine reduces the susceptibility of vaccinated individuals to infection, the reduction in infectiousness of vaccinated individuals and the reduction in infection fatality rate of vaccinated individuals. This new model can be described by the following system of equations:
We consider three cases. For comparison with the Base model, we consider the same level of efficacy: V Esym = V Esucpt = V Etrans = 65%. Next we consider a vaccine that is more effective at reducing the severity of illness but less effective at preventing infection and transmission. We consider a moderate case where V Esym = 70% and V Esucpt = V Etrans = 30%. Finally we consider a more extreme case where V Esym = 90% and V Esucpt = V Etrans = 10%. These three scenarios will be referred to as the “All 65%”, “70%-30%”, and “90%-10%” scenarios, respectively.
Iterations of main text Fig. 2 and Fig. 3 for these three alternative scenarios are shown further below. To compare and contrast cumulative vaccination results, directly below in Fig. 23 we show for the Base model and the three alternative scenarios the percentage of each group vaccinated after three months (left panel) and six months (right panel). In general, results are similar between the Base model with 65% of individuals 100% protected when vaccinated and the alternative where vaccines are 65% effective for all vaccinated. When shifts are apparent at either three or six months, they typically involve younger essential workers or younger seniors.
The cumulative percent of each demographic group (x-axis) vaccinated after three or six months under the Base model and three alternative scenarios for partial vaccine efficacy (y-axis) for each objective (panel). Vaccine supply is sufficient to vaccinate 10% of the population each month.
Across scenarios, we find that results after three months are sensitive in specific cases. When minimizing deaths or infections, results are relatively insensitive. A notable exception is for infections under the extreme scenario (90%-10%), we see a shift away from school-age children to the youngest and oldest. When minimizing YLL, we found substantial shifts between three groups: younger essential workers, younger seniors and older seniors.
After six months, the sensitivity and shifting lies mostly with younger groups (under 40) when focused on mortality (deaths or YLL) and conversely with older groups (over 40) when focused on infections (except for the extreme case, 90%-10%).
C.1 Vaccine reducing mortality and spread equally
C.2 Vaccine reducing mortality more than spread (moderately)
C.3 Vaccine reducing mortality more than spread (strongly)
Acknowledgements
This material is based upon work supported by the National Science Foundation (NSF) Graduate Research Fellowship Program under Grant No. 1650042. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF. MS is supported by an Emergency COVID-19 Seed Award from the California Breast Cancer Research Program of the University of California, Grant No. R00RG2419. GC is partially supported by NSF grant No.’s 2026797, 2034003, and NIH R01 GM 130900. We thank participants in the MIDAS network COVID-19 modeling seminar for helpful comments.
Appendix
Footnotes
↵1 We also note (13) use simulation without optimization to explore implications of vaccines with various levels of direct and indirect protection.
↵2 In the event that vaccination requires two doses over time, we consider an individual vaccinated upon receipt of the second dose at time t and we assume that v indicates the number of individuals that can be vaccinated with the required number of doses.
↵3 This vaccine efficacy is inclusive of any efficiency loss from typical handling in the distribution chain.
References
