Sequential allocation of vaccine to control an infectious disease.

The problem of optimally allocating a limited supply of vaccine to control a communicable disease has broad applications in public health and has received renewed attention during the COVID-19 pandemic. This allocation problem is highly complex and nonlinear. Decision makers need a practical, accurate, and interpretable method to guide vaccine allocation. In this paper we develop simple analytical conditions that can guide the allocation of vaccines over time. We consider four objectives: minimize new infections, minimize deaths, minimize life years lost, or minimize quality-adjusted life years lost due to death. We consider an SIR model with interacting population groups. We approximate the model using Taylor series expansions, and develop simple analytical conditions characterizing the optimal solution to the resulting problem for a single time period. We develop a solution approach in which we allocate vaccines using the analytical conditions in each time period based on the state of the epidemic at the start of the time period. We illustrate our method with an example of COVID-19 vaccination, calibrated to epidemic data from New York State. Using numerical simulations, we show that our method achieves near-optimal results over a wide range of vaccination scenarios. Our method provides a practical, intuitive, and accurate tool for decision makers as they allocate limited vaccines over time, and highlights the need for more interpretable models over complicated black box models to aid in decision making.

Introduction

The problem of optimally allocating a limited supply of vaccine to control a communicable disease has broad applications in public health: for example, in the control of diseases such as dengue, Ebola, seasonal influenza, and COVID-19. For many epidemic diseases such as COVID-19 and seasonal influenza, vaccination is an essential part of control; non-pharmaceutical and other interventions (e.g., treatment) may be helpful but insufficient. Often, however, vaccine supplies are constrained, so policy makers must address the question of how best to allocate limited vaccine supplies.

A number of studies have examined the general problem of optimal vaccine allocation. For the case of a one-time vaccine allocation, some studies use linear or mixed linear programming formulations with the objective of minimizing the number or cost of vaccines needed to reduce the epidemic reproduction number below 1 (e.g., [1], [2], [3]). Other studies develop simple analytical conditions to guide a one-time vaccine allocation, as a function of the objective to be optimized (e.g., minimize deaths, new infections, life years lost, quality-adjusted life years [QALYs] lost, or ) [4], [5]. For the case of vaccine allocation over time, some studies use optimal control formulations to assess the level of vaccination needed to reduce below 1 (e.g., [6], [7], [8]). One study uses a stochastic linear program to determine optimal allocation of limited vaccine doses to different geographic regions in two time periods [9]. A study of influenza vaccination showed that the optimal allocation is dynamic and depends on population structure: for example, in some cases it is optimal to first vaccinate high-transmission groups and then switch to vaccinating the most vulnerable groups in order to minimize deaths [10].

The vaccine allocation problem has received renewed attention during the COVID-19 pandemic [11]. For the case of a one-time allocation, some studies use numerical analysis of age-based and exposure risk-based compartmental models to evaluate alternative allocation policies (e.g., [12], [13], [14], [15]), generally finding that vaccinating older individuals minimizes deaths whereas vaccinating younger individuals or those in high-contact occupations minimizes transmission and new infections. Using numerical analysis of dynamic compartmental models, two studies show that dynamic allocation of COVID-19 vaccines to different population groups over time can be highly beneficial in minimizing infections, deaths, and life years lost, with the optimal allocation depending on the objective [16], [17]. For example, a study of COVID-19 in China showed that if the goal is to minimize new infections, then individuals aged 15–39 years old should be prioritized for vaccination until 47% of this group has been vaccinated, at which point it is optimal to prioritize individuals aged 40–64 until 26% coverage is reached; if the goal is to minimize deaths, 100% of individuals 65 and older should be vaccinated, followed by vaccination of individuals aged 40–65 until 97% coverage is reached [16].

As highlighted in existing studies, the vaccine allocation problem is highly complex and nonlinear, with no closed-form analytical solution, even in the simplest cases. Decision makers need a practical, accurate, and interpretable method to guide vaccine allocation. In this paper we develop simple analytical conditions that can guide the allocation of vaccines over time. A previous study that developed analytical conditions for one-time vaccine allocation among different population groups relied on a limiting assumption about the amount of vaccine that could be allocated to any group [4]. Here we extend that work to relax the upper bound on the amount of vaccine that can be allocated and then develop a method for allocating vaccines over time. The optimal allocation problem can no longer be reduced to a knapsack problem. Instead, we reduce the problem for any single time period to a piecewise linear optimization problem, and we then solve the multi-period problem considering the time periods sequentially. We illustrate our method with an example of COVID-19 vaccination, calibrated to COVID-19 spread in New York State.

Framework

Vaccine allocation problem

We consider an SIR model of a population with interacting groups in which an infectious disease is spreading. Individuals in each group can be susceptible (), infected (), recovered (), or dead (). Individuals in group can acquire infection from contact with individuals in their own population group (at rate ) or another population group (at rate ). Infected individuals in group either recover (at rate ) or die (at rate ). We consider a relatively short time horizon and thus do not include births, non-infection-related deaths, or other forms of entry into and exit from the population.

The compartmental model is governed by the following differential equations:

We assume that a preventive vaccine with effectiveness is available and that vaccination of susceptible individuals moves them to a recovered health state. Vaccination does not affect the transmission rates between infected and unvaccinated individuals () nor the recovery rates of infected individuals (). We let denote the (constant) population size.

For vaccine allocation, we consider a discrete set of time periods, and let , where denotes the amount of vaccine allocated to population group in time period . We assume that vaccination is instantaneous at the start of each time period. We let ; this is the set of vaccine allocations up through period . For notational simplicity we write the full set of vaccine allocations as . We assume that a limited number of vaccines, , are available to be distributed at the start of each time period .

We denote by , and the proportion of individuals in each compartment at time without vaccination. We let , , , and be the proportion of individuals in each compartment at time in the presence of vaccination where time falls within vaccination period . Let time denote the beginning of time period . We have The allocation instantaneously moves a fraction of population from the susceptible state to the recovered state at time . The epidemic evolves continuously and at discrete time points the trajectory is changed by vaccine allocation.

We consider four different objectives : cumulative new infections, deaths, life years lost, or quality-adjusted life years (QALYs) lost due to death, up to the end of time period . We write the vaccine allocation problem as follows: In each time period the available vaccine supply is limited to . Allocation to each population group must be nonnegative and is limited by the size of the susceptible (unvaccinated) population at the start of the time period.

Solution approach

Our goal is to determine the allocation of vaccines over time that minimizes . However, the vaccine allocation problem is highly nonlinear. Even during a single time period the epidemic evolves nonlinearly, allocations to one population group affect epidemic growth in the other population groups, and vaccine allocation can have a nonlinear effect on epidemic evolution. For this reason, we develop a solution approach in which we allocate vaccines in each time period based on epidemic conditions at the start of the time period.

For each time period , we thus solve the following optimization problem:

Single-period approximated optimal solution

A previous study that considered only a single time period used Taylor series approximations of the SIR model equations to develop optimality conditions for vaccine allocation for each of the four objectives we consider [4]. The problem reduces to a knapsack problem, where the coefficients of the knapsack problem are given by , , , and for the objectives of minimizing new infections, deaths, life years lost, and QALYs lost due to death, respectively. The term denotes the average number of life years lost for individuals in group who die from the disease and denotes the quality-of-life multiplier for individuals in group . The knapsack approach was shown to achieve optimal or near-optimal solutions for the case of a single time period [4]. However, the level of vaccination was restricted to be lower than , where is the length of the time period, and so that the approximations are always positive over the considered time horizon, and so that .

Here we develop optimality conditions for a single time period for higher levels of vaccination: in particular, we extend the analysis so that the upper bound on the level of vaccination is . The resulting problem is no longer a knapsack problem, but the objective functions are piecewise linear, which allows us to develop simple conditions characterizing the optimal solution. For simplicity in notation we omit the subscript and without loss of generality we let and denote the start and end of the time period, respectively.

For the objective of minimizing new infections, we approximate the disease dynamics at time using first-order Taylor expansions, and take the positive part so that the approximation of is positive for all :

For the objectives of minimizing deaths, life years lost, or QALYs lost due to death, we use second-order Taylor expansions of : Since we are considering a short period of time and do not consider entry into or exit from the population, is a trapping state. Therefore, the proportion of people who have died in each group has to be increasing over time. By construction of , the approximation of is indeed increasing for all .

Using these approximations, we derive the optimal vaccine allocation for each objective, assuming that . All proofs are in the Appendix. For notational simplicity, we will write . We define and .

Minimize New Infections. Once an individual has been infected, that person can either move to the recovered or dead compartment. Therefore, the proportion of the population that has been infected by time is given by The number of infections equals the sum over all population groups of the proportion of the population in the infected, recovered, and dead compartments (, , and , respectively) minus . The term corresponds to individuals who were effectively vaccinated (and thus transitioned to the recovered state) but who were never infected.

The approximated objective function is given by

This function is continuous and piecewise linear in .

If , the solution that minimizes satisfies . The term is the maximum level at which group can be vaccinated in an optimal solution.

Assume thatand that the population groups are ordered in decreasing order of their initial force of infection, i.e.  if . Define . The vaccine allocation that minimizes INF( In other words, it is optimal to vaccinate the groups in decreasing order of their initial force of infection.

Minimize Deaths. The approximated objective function is given by This function is a continuous piecewise linear function of .

If , the solution that minimizes satisfies . The term is the maximum level at which group can be vaccinated in an optimal solution.

Assume thatand that the population groups are ordered in decreasing order of their initial force of infection multiplied by the mortality rate, i.e.  if . Define . The vaccine allocation that minimizes is In other words, it is optimal to vaccinate the groups in decreasing order of their initial force of infection multiplied by the mortality rate.

Minimize Life Years Lost and QALYs Lost. We can write the functions for life years lost (LY()) and QALYs lost due to death (QALY()) as follows:

Since the functions and are weighted sums of , the solution to minimizing life years lost and QALYs lost due to death follows directly from the solution to minimizing deaths.

Assume that . To minimize , it is optimal to vaccinate groups in decreasing order of their initial force of infection multiplied by the mortality rate and expected life years lost ( ), up to the level . To minimize , the optimal solution is to vaccinate groups in decreasing order of their initial force of infection multiplied by the mortality rate and expected QALYs lost due to death ( ), up to the level .

For all four objective functions considered, the approximated optimal solution is an all-or-nothing allocation: we allocate as much of the vaccine as possible to groups in order of priority (up to the level for each group ) until no vaccine remains, and thus may allocate no vaccine to some population groups.

Dynamic allocation

We solve the dynamic vaccine allocation problem sequentially, applying the single-period optimality conditions to determine an allocation at the beginning of each time period based on the state of the epidemic at that point. The optimality conditions for all four objectives depend on the fraction of infected individuals when the allocation is made. We simulate the model numerically to determine for each time period . Algorithm 1 summarizes our solution approach.

Here we solve the optimal vaccine allocation problem sequentially using first- and second-order approximations. We note that even using first-order Taylor series expansions directly in the original multi-period problem (OPT1) yields a non-convex problem. For the case of a single time period, the resulting problem reduces to (OPT2), but for multiple time periods, the resulting problem is non-convex and can only be solved numerically; we provide an example for the case of two time periods in Supplemental Section B. We have instead chosen to use an approach that yields analytical insight.

Extensions

Fairness is often an important concern when allocating scarce medical resources such as vaccines [11], [30], [31], [32]. Our optimization problems (OPT1)–(OPT2) can be extended to include an equity constraint: where is the minimum fraction of the population in group that must be vaccinated in time period . We define the matrix such that . In this case, we consider as our decision variable with the constraints . With this equity constraint, the optimal solution follows the same priority order given by Propositions 2, 4, and 5, but each group is first allocated at least vaccines.

We can also modify the model to constrain the level of vaccination in group to be lower than , where is the vaccine acceptance rate in group . Such a constraint could reflect vaccine hesitancy or lack of vaccine access among certain populations. In this case, we define so that it is never optimal to vaccinate group above . Again, the optimal solution follows the priority order given by Propositions 2, 4, and 5, with each group being limited to vaccines in any time period.

In both cases our solution methodology still applies.

Example: COVID-19 vaccination

We illustrate our dynamic allocation method with the example of COVID-19 in New York State. We divide the population into age groups: individuals under age 20 (group 1), individuals aged 20–39 (group 2), individuals aged 40–65 (group 3), and individuals over age 65 (group 4).

Model instantiation and calibration

We instantiate the epidemic model using data for New York State from September 1, 2020 to November 30, 2020. We use daily COVID-19 cases and deaths in New York State, and obtain values for other model parameters from the published literature (Table 1).

Table 1

Values and sources for model parameters.

Parameter	Description	Value	Source
f_1	Fraction of individuals <20 years old	0.25	[18]
f_2	Fraction of individuals 20–39 years old	0.27	[18]
f_3	Fraction of individuals 40–65 years old	0.31	[18]
f_4	Fraction of individuals ≥65 years old	0.16	[18]
d_m	Average duration of mild infection (days)	11	[19], [20], [21], [22]
d_s	Average duration of severe infection (days)	8	[23], [24]
α_1	Fraction of infections that become severe for individuals <20 years old	0.02	[18], [25]
α_2	Fraction of infections that become severe for individuals 20–39 years old	0.15	[18], [25]
α_3	Fraction of infections that become severe for individuals 40–65 years old	0.26	[18], [25]
α_4	Fraction of infections that become severe for individuals ≥65 years old	0.46	[18], [25]
d_1	Average duration of infection for individuals <20 years old (days)	11.13	Calculateda
d_2	Average duration of infection for individuals 20–39 years old (days)	12.18	Calculateda
d_3	Average duration of infection for individuals 40–65 years old (days)	13.10	Calculateda
d_4	Average duration of infection for individuals ≥65 years old (days)	14.68	Calculateda
ξ_1	Infected fatality ratio for individuals <20 years old	0.0000988	[18], [26]
ξ_2	Infected fatality ratio for individuals 20–39 years old	0.0005750	[18], [26]
ξ_3	Infected fatality ratio for individuals 40–65 years old	0.0043939	[18], [26]
ξ_4	Infected fatality ratio for individuals ≥65 years old	0.0350831	[18], [26]
μ_1	Daily death rate for individuals <20 years old	0.0000088	Calculatedb
μ_2	Daily death rate for individuals 20–39 years old	0.000047	Calculatedb
μ_3	Daily death rate for individuals 40–65 years old	0.00034	Calculatedb
μ_4	Daily death rate for individuals ≥65 years old	0.00239	Calculatedb
γ_1	Daily rate at which individuals <20 years old recover and become immune	0.090	Calculatedc
γ_2	Daily rate at which individuals 20–39 years old recover and become immune	0.082	Calculatedc
γ_3	Daily rate at which individuals 40–65 years old recover and become immune	0.076	Calculatedc
γ_4	Daily rate at which individuals ≥65 years old recover and become immune	0.066	Calculatedc
η	Vaccine effectiveness	0.90	[27]
L_1	Expected life years lost for individuals <20 years old	69.29	[18], [28]
L_2	Expected life years lost for individuals 20–39 years old	50.28	[18], [28]
L_3	Expected life years lost for individuals 40–65 years old	29.81	[18], [28]
L_4	Expected life years lost for individuals ≥65 years old	12.95	[18], [28]
q_1L_1	Quality-adjusted expected life years lost for individuals <20 years old	63.02	[18], [28], [29]
q_2L_2	Quality-adjusted expected life years lost for individuals 20–39 years old	45.04	[18], [28], [29]
q_3L_3	Quality-adjusted expected life years lost for individuals 40–65 years old	27.50	[18], [28], [29]
q_4L_4	Quality-adjusted expected life years lost for individuals ≥65 years old	11.22	[18], [28], [29]

The average infection duration is calculated as .

The infected fatality ratio is estimated from the cumulative number of infections and deaths. Since we do not model the severity of the disease (mild vs. severe infection), the average death rate is calculated as .

The average recovery rate is calculated as .

We use model calibration to determine the transmission rate parameters () and the initial total number of infected individuals (). We assume that the initial distribution of cases is consistent with the age distribution; that is, for all . Since several studies have shown that the total number of cases could be many times higher than the number of confirmed cases [33], [34], we calibrate to a 7-day rolling average of reported deaths from September 1 to November 30, 2020 (Appendix Figure C.1a) and compare our model projections to multiples of a 7-day rolling average of new confirmed cases (Appendix Figure C.1b). We calibrate to daily deaths only up until December 1 since vaccination began on December 14 [35], and we want to capture the trend of the epidemic without vaccines.

We use Latin hypercube sampling for calibration, randomly sampling each parameter from a range of values [36]. We measure goodness of fit using the sum of squared errors. The calibrated transmission rate values are: and the starting initial infected population is . The resulting value is 1.14.

Fig. 1 compares the calibrated model’s output to the New York State data on deaths and confirmed cases. The model output closely matches the calibration target of reported deaths (Fig. 1(a)). The model’s projected total number of infected individuals is approximately 2.5 times higher than daily confirmed cases in New York State (Fig. 1(b)), which is consistent with studies such as [33], [34] that suggest that the total number of people infected could be 2–5 times the number of confirmed cases due to a large population of asymptomatic individuals and untested individuals.

Fig. 1

Calibrated model’s daily number of deaths, and multiples of daily confirmed cases compared to reported values (7-day rolling averages) for New York state.

As there is uncertainty about the number of COVID-19 cases, we also consider scenarios where there are two, five, and ten times [37] as many infected and recovered individuals as calibrated. We refer to the scenario that uses the calibrated values as the base case scenario.

Dynamic vaccine allocation

We consider three time periods of length: , or 30 days. We assume that a vaccine with effectiveness is available [27]. We assume that vaccines are available at the beginning of each time period ( and ). We consider available vaccination levels between 0% and 8% of the population in increments of 1%.

We have extended the analysis of a previous study [4] so that the upper bound on the level of vaccination is instead of ; in doing so, we now solve a piecewise linear optimization problem rather than a knapsack problem. Supplemental Table D.1 shows the maximum proportion of the susceptible population that can be vaccinated in the first time period for each scenario and time horizon considered. For short time horizons ( days), the bound of corresponds to more than 94.8% of the susceptible population, compared to only 16.5% with the previous bound of . For longer time horizons ( days), the bound of corresponds to approximately 55% of the susceptible population, compared to only 8.2% with the previous bound of .

Using the parameter values from the calibration and the literature (as detailed in Section 3.1), we determine the optimal vaccine allocation using Propositions 2, 4, and 5 for each time horizon, time period, vaccination level, and objective function. Table 2 summarizes our findings. The optimal vaccine allocation did not vary across scenarios (when varying the initial proportion of infected and recovered individuals) so we only present results for the base case scenario.

Table 2

Optimal vaccine allocation for each time horizon, time period, vaccine level, and objective function, as determined using the dynamic allocation method.

Length of time periods	Time period	Available vaccine level	Optimal group order of vaccination
			Minimize new infections	Minimize deaths, LYs and QALYs lost
T=7 days	τ=1	0%≤NP≤8%	2, 1, 3, 4	4, 3, 2, 1
	τ=2	0%≤NP≤8%	2, 1, 3, 4	4, 3, 2, 1
	τ=3	0%≤NP≤6%	2, 1, 3, 4	4, 3, 2, 1
		7%≤NP≤8%	1, 2, 3, 4	4, 3, 2, 1
T=15 days	τ=1	0%≤NP≤8%	2, 1, 3, 4	4, 3, 2, 1
	τ=2	0%≤NP≤8%	2, 1, 3, 4	4, 3, 2, 1
	τ=3	0%≤NP≤5%	2, 1, 3, 4	4, 3, 2, 1
		6%≤NP≤7%	1, 2, 3, 4	4, 3, 2, 1
		NP=8%	1, 3, 2, 4	4, 3, 2, 1
T=30 days	τ=1	0%≤NP≤8%	2, 1, 3, 4	4, 3, 2, 1
	τ=2	0%≤NP≤8%	2, 1, 3, 4	4, 3, 2, 1
	τ=3	0%≤NP≤4%	2, 1, 3, 4	4, 3, 2, 1
		NP=5%	1, 2, 3, 4	4, 3, 2, 1
		6%≤NP≤8%	1, 3, 2, 4	4, 3, 2, 1

For the objective of minimizing infections, group 2 is vaccinated first for any time period length and available vaccine level in the first and second time periods. The allocation switches to a different group in the third time period if the available vaccine level is high enough: when the length of the time periods is short ( days), the allocation switches to vaccinating group 1 before group 2 if ; when the length of the time periods is longer ( days), group 3 is vaccinated before group 2 and after group 1 for high enough levels of vaccine (for when days, and for when days).

To minimize deaths, life years lost, and QALYs lost due to death, the order of vaccination is group 4, followed by group 3, then group 2, and finally group 1. This result holds across all time periods, time period lengths, and available vaccine levels. The result is intuitive because mortality increases significantly with age () and the expected gain in life years or QALYs from vaccinating younger age groups is not enough to offset the higher mortality of older age groups.

Quality of decisions

To evaluate the accuracy of the approximated optimal allocations, we compare the solutions to allocations determined using the exact Eqs. (1). We determine the true optimal solution via exhaustive search, discretizing the range of feasible vaccine allocations ( and ) using a grid size of 0.001. We present results for each time period length: , or days.

7 days.Fig. 2 compares the value of the objective function at the numerical and approximated optimal solution when days.

Fig. 2

days. Percentage difference between the approximated and numerical optimal solutions.

For the objectives of minimizing deaths, life years lost, and QALYs lost due to death, the approximated optimal solution matches the numerical optimal solution found via exhaustive search in all cases (for any level of vaccination and time period).

For the objective of minimizing new infections, the approximated and numerical solutions differ slightly in the third time period (for ), but the percentage difference in new infections is small (). The discrepancy in allocations occurs close to the switching point (Fig. 3): the approximated solution allocates all vaccines to group 2 when , and to group 1 when . In contrast, the numerically optimal solution allocates vaccines to both groups 1 and 2 when . The discrepancy occurs because the approximated solution method leads to an all-or-nothing allocation, which is not always optimal.

Fig. 3

days, . Numerical and approximated optimal vaccine allocation to minimize new infections.

15 days.Fig. 4 compares the value of the objective functions at the numerical and approximated optimal solution for a time horizon of days.

Fig. 4

days. Percentage difference between the approximated and numerical optimal solutions.

For the objectives of minimizing deaths, life years lost, and QALYs lost, the approximated and numerical optimal solution match in all cases.

For the objective of minimizing new infections, the approximated and numerical solutions differ in the second time period when , and in the third time period when . However, the percentage difference in new infections is small (). The discrepancy between the two solutions occurs near a switching point (Fig. 5). For example, in the third time period, the approximated solution allocates vaccines to group 2 followed by group 1 when , to group 1 followed by group 2 and then group 3 when , and to group 1 followed by group 3 when . The numerical optimal solution allocates vaccines to both groups 1 and 2 when , to groups 1, 2 and then 3 when , and to groups 1 and 3 when . The switch occurs more gradually in the numerical solution because the approximated solution is an all-or-nothing allocation.

Fig. 5

days. Numerical and approximated optimal vaccine allocation to minimize new infections.

30 days. Finally, Fig. 6 compares the value of the objective functions at the numerical and approximated solution for days. For deaths, life years, and QALYs, the objective function values vary by at most 0.005%, and for new infections, the objective function values vary by at most 0.02%. The solutions match for lower levels of vaccination (, , in the first, second, and third time periods, respectively) but diverge for higher levels of vaccination. Again, this is because the approximate solution in each period is an all-or-nothing allocation whereas the numerical optimal solution involves a partial allocation between groups (Fig. 7). The allocations also differ due to upper limits on allocation in the approximate method: for example, in time period for , when considering the objective of minimizing deaths, life years lost, or QALYs lost, the numerical optimal solution allocates all vaccines to group 4 whereas the approximate solution allocates vaccines to both group 4 and group 3.

Fig. 6

days. Percentage difference between the approximated and numerical optimal solutions.

Fig. 7

days. Numerical and approximated optimal allocations.

Supplemental Table D.2 summarizes the percentage difference between the approximated and numerical optimal solution for each time horizon, time period, vaccine level, and objective function. The approximated and numerical solution differ by at most 0.018%. Supplemental Table D.3 shows the run times in seconds to find the approximated and numerical optimal solution for all vaccination levels and objective functions, for each scenario, time horizon, and time period. The average run time for the approximated solution was no more than a tenth of a second whereas the average run time for the numerical solution was on the order of hundreds or thousands of seconds.

Discussion

Our method for dynamic allocation of limited vaccines to control an infectious disease uses a simple rule based on epidemic conditions at the start of each time period. Using numerical simulations of COVID-19 in New York State we show that our method achieves near-optimal results over a wide range of vaccination scenarios.

The approximated optimal solution is an all-or-nothing allocation. Numerical simulations show that the approximated and numerical optimal solutions match for lower levels of vaccination. For higher levels of vaccination where partial allocation is optimal, we find that discrepancies occur near switching points, but the numerical solution nevertheless follows the priority order given by the approximated solution, and the numerical and approximated optimal solution differ by at most 0.018%. In practice, when considering longer time horizons, the accuracy of the allocation method can be improved by considering more time periods with shorter lengths.

The method we developed is practical for decision makers as it provides simple, interpretable conditions for vaccine allocation over time, and avoids the need to solve an intractable nonlinear dynamic optimization problem. Moreover, the simple analytical conditions match our intuition: to minimize new infections in any period it is optimal to vaccinate groups in decreasing order of their force of infection, whereas to minimize deaths it is optimal to vaccinate groups in decreasing order of their force of infection multiplied by the mortality rate. For the example of COVID-19, we find that individuals aged 20–39 should be prioritized for vaccination in the first time period to minimize new infections because they have the highest initial force of infection. However, prioritizing older individuals (over age 65) minimizes deaths, life years lost, and QALYs lost due to death because of the high mortality rate in this group.

Our method is also practical because it generates a prioritized list of population groups for vaccine allocation: for instance, in the COVID-19 example, to minimize new infections in the first time period, the priority order for vaccination is first individuals aged 20–39, then individuals under age 20, then individuals aged 40–65, and finally individuals over age 65. In practice, public health decision makers frequently allocate resources using prioritization schemes of this type. For example, several U.S. states adopted age-based COVID-19 vaccine allocation guidelines [38], [39], [40], [41]. As another example, the Centers for Disease Control and Prevention’s (CDC’s) recommendations for COVID-19 vaccine allocation prioritize groups based on age, underlying medical conditions, and occupation [42].

Our analysis has several limitations. We assumed that all available vaccine doses in each time period are used. However, the number of individuals receiving vaccination may vary due to factors such as access to vaccination or vaccine hesitancy [11], [43]. Additionally, epidemic growth is typically stochastic while our method uses a deterministic model to project epidemic evolution. Since our method determines vaccine allocation at the beginning of each time period, the model could be updated over time as more information becomes available. For example, instead of forecasting the fraction of infected people at the start of each time period using numerical simulations, policy makers could observe actual infection trends and estimate the fraction of infected people in each population subgroup. Additionally, the model parameters could be recalibrated at the start of each time period before generating the resulting allocation. Further work could automate our decision tool and could extend the analysis to a stochastic setting.

The problem of allocating limited vaccines is important, but complex. Our allocation method provides a practical, intuitive guide for decision makers that can achieve near-optimal solutions as they allocate limited vaccines over time. For the COVID-19 example, our findings are consistent with CDC and other recommendations to prioritize older age groups for vaccination so as to minimize deaths [11], [42]. For vaccination against other communicable diseases, our method can also provide practical, intuitive solutions. Although black box models are prevalent in the literature on vaccine allocation, our study shows that accuracy need not be sacrificed for interpretability. Our analysis highlights the need for interpretable models to aid in important problems in public health and epidemic control, including that of allocating limited vaccines to control a communicable disease.

Declaration of Competing Interest

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