Abstract
A vaccine, when available, will likely become our best tool to control the current COVID-19 pandemic. Even in the most optimistic scenarios, vaccine shortages will likely occur. Using an age-stratified mathematical model, we determined optimal vaccine allocation for four different metrics (deaths, symptomatic infections, and maximum non-ICU and ICU hospitalizations) under a wide variety of assumptions. We find that a vaccine with effectiveness ≥50% would be enough to substantially mitigate the ongoing pandemic provided that a high percentage of the population is optimally vaccinated. When minimizing deaths, we find that for low vaccine effectiveness, it is optimal to allocate vaccine to high-risk (older) age-groups first. In contrast, for higher vaccine effectiveness, there is a switch to allocate vaccine to high-transmission (younger) age-groups first for high vaccination coverage. While there are other societal and ethical considerations, this work can provide an evidence-based rationale for vaccine prioritization.
Introduction
As of 13 August 2020, over 750 thousand people have died due to the ongoing SARS-CoV-2 pandemic (1). Different countries have enacted different containment and mitigation strategies, but the world awaits impatiently for the arrival of a vaccine as the ultimate tool to fight this disease and to allow us to resume our normal activities. There are over 100 vaccines under development (2,3), with some currently undergoing phase 3 clinical trials (3). However, there are many unknowns surrounding a potential vaccine, including how efficacious it would be, how long it would be protective, how effective it would be in older individuals, how many doses would be immediately available and how long scaling up the vaccine production would take. Furthermore, should early vaccines have low effectiveness, what are the potential trade-offs between using a low-effectiveness vaccine and waiting for a vaccine with a more desirable vaccine effectiveness? With the hope of producing a vaccine in the near future comes the difficult task of deciding who to vaccinate first as vaccine shortages are inevitable (4–6). Here we utilized a mathematical model paired with optimization algorithms to determine the optimal use of vaccine for 100 combinations of vaccine effectiveness (VE) and number of doses available under a wide variety of scenarios.
Results
Briefly, we developed a deterministic age-structured mathematical model of SARS-CoV-2 transmission with a population stratified into 16 age-groups (Fig. S1, SM). Because, historically, vaccine is distributed to each state in the United States proportional to its population, and the allocation strategy is then determined at the state level (7), we chose a state level model with a population similar to Washington State in size and demographics; however, our results are generalizable to other populations. We assumed that children were less susceptible to infection than middle-aged adults (20 to 65 years old), while older adults (older than 65) were relatively more susceptible (8). We assumed that both natural and vaccine-induced immunity last at least one year (our time horizon). At the beginning of our simulations, 20% of the population have already been infected and are immune (additional results for 10, 30 and 40% of the population can be found in the SM) and all social distancing interventions have been lifted. Here, we consider that front-line health care workers, who should obviously be prioritized, have already been vaccinated.
For the vaccine optimization, we collated the 16 age-groups into five vaccination groups: children (aged 0–19), adults between 20 and 49 years old, adults between 50 and 64 years old, adults between 65 and 74 years old, and those 75 and older. This stratification reflects our current knowledge of disease severity and mortality based on age (9, 10). We developed an optimization routine that combined a coarse global search algorithm with a fast optimizer to explore the entire space of possible combinations of vaccine allocation. We compared the optimal allocation strategy given by the optimizer to a pro-rata allocation, where the vaccination coverage to each vaccination group is distributed proportionally to population size in each group. We considered VE ranging from 10% to 100% and vaccination coverage ranging from 10% to 100% of the total population. We evaluated four objective functions reflecting different metrics of disease burden that could be considered by decision makers: minimization of the total number of symptomatic infections, total number of deaths, number of cases requiring hospitalization (non-ICU) at the epidemic peak, and number of cases requiring ICU hospitalization at the epidemic peak. The last two objective functions were chosen because hospital bed (non-ICU and ICU) occupancy is a key metric currently used to determine county/state/country readiness to move between different interventions strategies. Here, we utilized the total number of licensed ICU beds in WA state and its current goal of staying below 10% of hospital beds occupied by COVID-19 cases (11, 12) as references when interpreting our results.
Epidemic mitigation and containment: Our model suggests that herd immunity will be achieved once 60% of the population is infected (equivalently 40% vaccinated with a perfect vaccine assuming 20% of the population has already immunity) (Fig. 1J, Fig. S2 and Fig. 2A).
The epidemic can be substantially slowed with any vaccine with a VE ≥ 50% as long as a majority of the population is vaccinated (Fig. 1E, Fig. 2A), and over 50% of deaths could be averted with as little as 35% of the population vaccinated (Fig. 2A, B). If VE = 60%, the epidemic is completely contained if we optimally vaccinate 70% of the population, but we would only need to vaccinate 50% of the population if VE = 70% (Fig. 2A and Fig. 1F, I). Only vaccines with VE ≥50% can maintain the number of non-ICU hospitalizations below the established goal (≤10% hospital-bed occupancy by COVID-19 patients) and can prevent an overflow of the ICUs. With VE = 60%, 54% of the population would have to be optimally vaccinated to satisfy both conditions (Figs. 2C, D, S3F and S4F). In contrast, with the same VE, over 67% of the population would have to be vaccinated with the pro-rata allocation in order to maintain hospitalizations (both non-ICU and ICU) below the desired goals (Figs. S6, S7, S8, S9). Utilizing the optimal allocation strategy matters most when less vaccine is available, with a maximum difference of 32% deaths averted (for VE = 100% and with enough vaccine to cover 20% of the population) and 32% symptomatic infections averted (for VE = 60% and vaccination coverage of 60%) when compared with a pro-rata allocation strategy (Fig. 3, S5, S10). As VE increases, both strategies tend to perform similarly as vaccination coverage increases (Fig. 3, S5 S10).
Optimal vaccine allocation changes with VE and vaccination coverage
The optimal allocation strategy to minimize total deaths is identical for vaccine efficacies between 10% and 50%: with low vaccination coverage, it is optimal to allocate vaccine first to the highest risk group (people over 75 years old) and then to the younger vaccination groups as more vaccine becomes available (Fig. 4A–E). However, there is a threshold phenomenon observed for VE≥60%: for low coverage, the optimal allocation is still to vaccinate the high-risk groups first, but when there is enough vaccine to cover roughly half of the population (60% for VE = 60%, 50% for VE = 70% and 40% for VE≥ 80%), there is a switch to allocate vaccine to the high-transmission groups (those aged 20–50 and children) first. This is because directly vaccinating those who are driving the epidemic results in a much slower epidemic curve and hence in fewer deaths (Fig. 1F–H). As more vaccine becomes available the optimizer allocates it to high-risk groups again (Fig. 4F–J).
Optimal vaccine allocation differs for different objective functions
Next, we investigated how the optimal allocation strategy changed for different objective functions and present results for VE = 60%. The optimal vaccine allocation for the four objectives were most different when fewer vaccines are available (enough vaccine to cover less than 30% of the total population). When minimizing symptomatic infections and peak non-ICU hospitalizations, priority was given to the younger vaccination groups, as they have the most contacts in our model and hence drive transmission (Fig. S11A,B). As we move toward more severe outcomes (ICU hospitalizations at peak and deaths), for which older individuals are most at risk, the optimal allocation strategy shifts toward those vaccination groups (Fig. S11C,D). Once more vaccine becomes available, the optimal allocation strategies are very similar for all objective functions. In fact, they are nearly identical for all the objective functions when there is enough vaccine to cover 60% and 70% of the population. For high coverage, the optimal allocation strategies for all objective functions shifted towards the high-transmission groups. To note, we did not impose the optimizer to use all the available vaccine. As a result, the optimizer found allocation strategies utilizing less than the total vaccine available while performing equally well. This was very prominent when VE and vaccination coverage were very high. For example, when minimizing peak ICU hospitalizations and VE = 90%, the optimizer utilized vaccine to cover 75% of the population even though there was vaccine available to cover the entire population. This is expected, complete containment is attained once a high proportion of the population is vaccinated and any vaccine used above that threshold will result in the same mathematical outcome.
Optimal vaccine allocation as a function of pre-existing immunity
As the COVID-19 pandemic dynamics have been dramatically different in different parts of the world, we expect to see a range of population-level naturally-acquired immunity when vaccination campaigns start. Hence, we investigated the optimal use of vaccine with 10%, 30% and 40% of the population is already immune at the beginning of the simulations. For all of these, the same pattern is observed when minimizing deaths: for low coverage, it is optimal to allocate all of the vaccine to the high-risk groups, for higher coverage, the optimal vaccination strategy switches to allocate more vaccine to the high-transmission groups. This threshold however varies with the degrees of pre-existing immunity in the population. When only 10% of the population is immune, the switch occurs when 80% of the population is vaccinated, but this threshold is observed when 40% of the population is vaccinated if 40% of the population has infection-acquired immunity prior to vaccination (Fig. S12). In addition, under low pre-existing immunity, the optimal strategy favors more the older vaccination groups (Fig. S12A), while under higher pre-existing immunity, the optimal allocation strategy tends to distribute vaccine more evenly across vaccination groups (Fig. S12D).
Robustness of optimal allocation strategies around major parameters: We explored the robustness of the optimal allocation strategies around key parameters of the transmission and natural history of SARS-CoV-2. Because susceptibility to SARS-CoV-2 infection remains un-clear, we compared the optimal allocation strategy under the assumption of differential susceptibility, as suggested in (8, 13) (presented throughout the text), to one assuming equal susceptibility across age-groups (Figs. S13, S14, S15 and S16), as suggested in (14, 15). The optimal allocation strategies under both equal and differential susceptibility were remarkably consistent for minimizing deaths, with more vaccine allocated to children as VE increases (VE≥60) and more vaccine becomes available (coverage to vaccinate 70% or higher) (Figs. 4 and S13). When minimizing symptomatic infections, optimal allocation strategies under equal susceptibility tended to allocate more vaccine to children than did the strategies under differential susceptibility (Figs. S17 and S14). The major differences were observed when minimizing peak hospitalizations (both non-ICU and ICU). Assuming equal susceptibility resulted in optimal allocation strategies that favored the high-transmission groups (children and adults aged 20 to 49 years old) over the high-risk groups (Figs. S18, S15, S19 and S16). In addition, we selected four parameters for which there is the most uncertainty and re-ran the optimization routine for several combinations (full details in SM). The optimal allocation strategies were very robust under this analysis (Supplemental Files. SF1–SF4). Finally, we compared the optimal allocation strategies when the simulations were started with a higher number of infected individuals (10,000 current infections). This would reflect a situation where the epidemic is in full exponential growth when vaccination becomes available. The optimal allocation strategy was surprisingly robust under this scenario, with nearly identical allocation strategies for all objective functions (Figs. S20, S21, S22, S23).
Discussion
The COVID-19 pandemic has devastated families and societies around the world. A vaccine, when available, would most likely become our best tool to control the spread of SARS-CoV-2. However, in the short term, even in the most optimistic scenarios, vaccine production would likely be insufficient. In this work, we paired a mathematical model of SARS-CoV-2 transmission with optimization algorithms to determine optimal vaccine allocation strategies. Given the current uncertainties surrounding such a vaccine (we do not know yet if and when this vaccine would be available, how efficacious it will be and the number of doses immediately available) we explored 100 combinations of VE and vaccination coverage under a wide variety of scenarios minimizing four metrics of disease burden.
Our results suggest that any vaccine with medium to high effectiveness (VE ≥ 50%) would be able to considerably slow the epidemic while keeping the burden on healthcare systems manageable, as long as a high proportion of the population is optimally vaccinated. Moreover, once VE = 70%, full containment of the epidemic would be possible. This is in agreement with vaccine modeling studies (16,17). Further, we showed that much can be achieved even with low vaccination coverage; indeed, with medium VE, over half of deaths can be averted by optimally vaccinating only 35% of the population. When minimizing deaths, for low VE and a low supply of vaccine, our results suggest that vaccines should be given to the high-risk groups first. For high VE and high vaccination coverage, the optimal allocation strategy switched to vaccinating the high-transmission groups (younger adults and children). This remained true under equal or reduced susceptibility to infection for children, pointing to the importance of children as key players in disease transmission. This finding is consistent with others (18, 19) that have found for other respiratory viruses that protecting the high-transmission groups can be the optimal use of resources.
Here, we utilized mathematical optimization to determine the optimal vaccine allocation and by design, did not impose any restrictions in the allocation strategies. However in practice, implementation of optimal strategies must also account for other factors (ethical, political and societal). When large quantities of vaccine are available, a feasible solution could involve first vaccinating the high-risk groups and then allocating the remaining vaccine to the high-transmission groups.
This study has several limitations. Our model assumes that both naturally- and vaccine-acquired immunity will last for at least one year. We do not yet know how long immunity against SARS-CoV-2 will last and there is some evidence that neutralizing antibodies become undetectable after just a few weeks following infection (20), though it is unclear how this correlates with immunity. If immunity were short-lived, then these results would only be applicable for that duration. Further, we assumed that asymptomatic and symptomatic infections would confer equal immunity. However, some studies have suggested that asymptomatic infections might result in a weaker immune response (21). We utilized mortality and hospitalization rates that were based on the epidemic in Wuhan, but these rates may vary vastly in different regions. Further, we compared modeled peak hospitalizations to current state goals for hospital bed occupancy, but deterministic models tend to overestimate the transmission dynamics, so it is possible that a lower vaccine effectiveness or a lower vaccination coverage could achieve the same goals. To keep the optimization from being unreasonably long, our model does not capture geographical differences or other heterogeneities. We assumed a vaccine that would only reduce susceptibility to infection, but other effects, e.g., a reduction in disease severity, might occur. We have identified optimal allocation strategies and once more information about a vaccine characteristics is known, validating our allocation strategies with more complex models is welcome. To avoid confounding effects from different interventions, we optimized vaccine allocation assuming no social distancing interventions in place. In reality, vaccination, at least at the beginning, would take place while some social distancing interventions remain in effect. Under those circumstances, we would need less vaccine to control the epidemic. In that sense, our results are conservative. We computed the optimal allocation strategies utilizing age as the sole risk factor. However, several studies (22) have shown that, as a result of health systems with systemic health and social inequalities, people from racial and ethnic minority groups are at increased risk for getting sick and dying from COVID-19 in certain countries. This is a crucial consideration that will be included in further studies and can point towards who, within a given age-group, should get the vaccine first.
We believe that these results can provide a quantification of the effectiveness of different allocation scenarios under four metrics of disease burden and can be used as an evidence-based guidance to vaccine prioritization.
Data Availability
Data is available in the manuscript.
Funding
LM, TL and ERB acknowledge support from NIAID, grant 1 UM1 AI148684–01.
Author contributions
LM and ERB designed the research. LM, and JE performed the research. LM and TL analyzed the data. LM, JE, TL and ERB wrote the manuscript.
Competing interests
LM and TL have received funding from Wellcome Trust for research unrelated to this manuscript. All other authors declare no competing interests.
Data and materials availability
Data are available in the main text or supplementary materials. Code available at https://github.com/lulelita/vaccine_optimization.
Supplementary materials
Materials and Methods
Mathematical model
We developed a deterministic mathematical model with 16 age-groups: 0–4, 5–9, 10–14, 15– 19, 20–24, 25–29, 30–34, 35–39, 40–44, 45–49, 50–54, 55–59, 60–64, 65–69, 70–74, and ≥75. For each age-group, the population is divided into the following compartments: those who are susceptible (S) to infection; exposed (E) but are not yet infectious; infectious (I); and recovered (R). Infectious individuals are classified by severity of symptoms. Those who are infectious may be asymptomatic (IA) or pre-symptomatic (P). As pre-symptomatic individuals become symptomatic, they may not require hospitalization (IS), require hospitalization (IH), or require intensive care (IC). Because of the short duration of our simulation, we did not model any births or deaths. We assumed a population of 7.615 million people, the current population of Washington State (23) and parameterized the demographic composition of the population to be similar to that of the United States (23).
The model equations (where a dot represents differentiation with respect to time) are: where the force of infection is
A diagram of the model is shown in Fig. S1.
Transmission between a susceptible and an infectious individual occur using an age-specific contact matrix (24) at an assumed R0 = 3. We used previously reported age-specific estimates of the severity of symptomatic cases that require hospitalization, critical care, and that lead to death (10). An average hospital stay lasts 5 days without intensive care, or 10 days if intensive care is required (9). We assumed a total duration of infectiousness of non-hospitalized infections of 5 days (14, 25). We assumed that immunity would last at least one year, which corresponds to the duration of our simulations.
A proportion, (1 − a), of infections is symptomatic. Of these, a proportion, h, requires hospitalization. Further, a proportion, c, of hospitalizations requires intensive care (ICU). Symptomatic and asymptomatic individuals are assumed equally infectious. Pre-symptomatic individuals, however, are assumed to be more infectious to match the finding that 40% of transmission occurs prior to symptom onset (9, 26). Once hospitalized, individuals are assumed no longer infectious. We determined the optimal allocation strategies assuming both that susceptibility to infection is age-specific and that susceptibility to infection is equal across all age-groups. See Table S1 for details.
Simulations were run with initial conditions set to a 20% recovered population and 1,000 current infections, distributed proportionally to population size (pro-rata) and disease severity, respectively.
Vaccination
We assumed a leaky vaccine (27) that reduces susceptibility to infection, modeled as a reduced probability of acquiring a SARS-CoV-2 infection. For each vaccination coverage and vaccination strategy considered, we computed within each age-group the fraction of susceptible people among all those individuals in that group who could have sought the vaccine (susceptible, infected pre-symptomatic, infected asymptomatic, and recovered asymptomatic populations), and utilized that fraction as the fraction of people who were actually vaccinated in each age-group, while assuming that the remaining vaccine would be wasted.
Uncertainty analysis
We performed uncertainty analysis in two ways. First, we examined the robustness of the optimal vaccination strategies against the number of initial infections, the fraction of the population with natural immunity at the start of the simulations and the relative susceptibility of children and older adults to infection. In addition, we selected four parameters (the basic reproduction number R0, the fraction of asymptomatic infections and the reduction/increase in infectiousness for asymptomatic and pre-symptomatic infections) that we considered would have the most influence in the optimal allocation strategies and re-run the optimization for the 100 VE-vaccination coverage pairs with 36 different combinations of those parameters. Table S2 provides the combination of values used for this part of the analysis.
Second, we examined the uncertainty in the output measures (number of infections, number of deaths, etc) arising from uncertainty surrounding the model parameters, and chose those parameters for which there is the least agreement. These parameters were: R0, the mean presymptomatic period, the mean infectious period of non-hospitalized symptomatic infections, the proportion of infections that are asymptomatic, and the relative infectiousness of asymptomatic infections. We sampled 1000 parameter sets from pre-determined distributions as follows. For R0, the proportion of asymptomatic infections, and the multiplier increasing or decreasing the infectiousness of pre-symptomatic infected individuals, we used truncated normals. For the duration of the latent and infectious periods (both pre- and post symptoms) we utilized gamma distributions. Finally, for the relative infectiousness of asymptomatic infections, we chose among three values: 0.2, 0.5 and 1. Table S1 gives the ranges utilized for each of these parameters. We sampled 1000 parameter sets and ran the model with those sets for the optimal allocation strategy for each combination of VE and vaccination coverage pair. Further, we removed the top and bottom 2.5% of the simulations.
Optimization
Objective functions
We performed the optimization routine to minimize four different objective functions: total symptomatic infections, total deaths, maximum number of hospitalizations not requiring intensive care and maximum number of hospitalizations requiring intensive care. For each of these, we ran the deterministic model for one year. Further, we evaluated the optimal vaccine allocation for 100 VE-vaccination coverage combinations (VE ranging from 10 to 100% and vaccination coverage ranging from 10 to 100% of the total population, in increments of 10%).
Feasibility
For the optimization routine, we collapsed the population into 5 vaccination groups: 0–19 years old, 20–49 years old, 50–64 years old, 65–74 years old and those aged 75 and older. We defined a decision variable in terms of the proportion total vaccine available: let v be a decision vector v = (v1, v2, … , v5), where vi represents the fraction of the total available vaccine to be given to vaccination group i. These vaccination groups were chosen for two reasons; clinically, these groups represent the finest granularity that we could find in the data and computationally, this reduces the dimension of the optimization problem.
Defining the decision variable in terms of proportion of total vaccine available allowed us to perform an exhaustive search on a coarse grid of the decision variable space that is independent of the population in each group. Knowing the total number of vaccines available let us seamlessly translate the decision variable between the fractions of vaccine to be used in each vaccination group and fractions of the population in each vaccination group to be vaccinated. In this setup, a feasible vector is defined as any vector v so that . This guarantees that we are using at most all of the available vaccine. However, this will not guarantee that the amount of vaccine given in each group does not exceed the actual population in that group. Hence, for any given total vaccine available, we further “repaired” any given feasible vector by mapping the feasible vector v into a feasible vector x = (x1, x2, … , x5) where xi represents the fraction of the population in vaccination group i to be vaccinated and removing any excess of vaccine for each group (so that for each xi, 0 ≤ xi ≤ 1). In this way, we can guarantee feasibility in terms of both the percentage of vaccine available and people in each vaccination group. This approach to modeling the decision variable and repairing a vector if needed simplified our application of an unconstrained optimization method (28) to this constrained problem by making the feasible set the unit simplex in dimension 5.
To run the deterministic model described above, the decision variable is then transformed into a vector of size 16, where each entry corresponds to one of the 16 age-groups of our model. Vaccine in each age-group was distributed proportionally to the fraction of the population that it represents within its vaccination group. For example, in the first group (people under 20 years old), children aged 0–4, 5–9, 10–14 and 15–19 represent 24%, 24.4%, 25.7%, and 25.9% of that group. Hence, within the first vaccination group, vaccine would be allocated to each of the four age-groups according to those proportions. The objective functions were then evaluated via the mathematical model.
Optimization routine
Our optimization routine consisted of two steps: first, we performed an exhaustive search on a coarse grid (29) of the unit simplex in the vaccination group space (the set of vectors (v1, v2, … , v5) with non-negative entries such that ). The grid was chosen so that the unit simplex was divided into 0.5 units and was computed in Sage (30). For example, a point in this grid is v = (0.05, 0.5, 0.1, 0.2, 0.15), that corresponds to utilize 5, 50, 10, 20 and 15% of the available vaccine in vaccination groups 1 through 5 respectively (as noted above, this vector would be repaired if vaccine exceeds the population in any vaccination group). For each point in the coarse grid, the four objective functions were evaluated. For each of the four objective functions, we selected the best 25 decision variables obtained in the grid search, the pro-rata allocation vector and an additional 25 decision variables sampled uniformly from the unit simplex (31), and used these 51 points as initial points for the Nelder-Mead minimizer implemented in SciPy (28, 32). To note, we also ran the optimization using particle swarm optimization algorithm (33, 34) and obtained similar optimal solutions, but solutions obtained with Nelder-Mead gave slightly better optimal values and ran faster.
Sampling from the entire simplex enabled exploration of the whole decision variable space to prevent the optimizer from getting stuck at a local minimum, which is common when utilizing this type of heuristic algorithm. The mathematical model and the optimization routine were implemented in Python.
Supplemental Files captions
SF1: Robustness analysis for the optimization strategies when minimizing symptomatic infections. Each page of the document represents one VE. Within each page, each subplot represents the optimal allocation strategy for that particular VE and for the combination of parameters given in the row of Table S2 indicated in the subtitle.
SF2: Robustness analysis for the optimization strategies when minimizing total deaths. Each page of the document represents one VE. Within each page, each subplot represents the optimal allocation strategy for that particular VE and for the combination of parameters given in the row of Table S2 indicated in the subtitle.
SF3: Robustness analysis for the optimization strategies when minimizing peak non-ICU hospitalizations. Each page of the document represents one VE. Within each page, each subplot represents the optimal allocation strategy for that particular VE and for the combination of parameters given in the row of Table S2 indicated in the subtitle.
SF4: Robustness analysis for the optimization strategies when minimizing peak ICU hospitalizations. Each page of the document represents one VE. Within each page, each subplot represents the optimal allocation strategy for that particular VE and for the combination of parameters given in the row of Table S2 indicated in the subtitle.
Supplemental Figures
Supplemental Tables
Acknowledgments
LM acknowledges Mia Moore for helpful discussions regarding the model structure and Michael Gutteridge for help with cluster computing.