Where to locate COVID-19 mass vaccination facilities?
Dimitris Bertsimas, Vassilis Digalakis Jr., Alexander Jacquillat, Michael Lingzhi Li, Alessandro Previero
WWhere to locate COVID-19 mass vaccination facilities?
D. Bertsimas, V. Digalakis Jr., A. Jacquillat, M. L. Li, A. Previero
Sloan School of Management and Operations Research Center, Massachusetts Institute of Technology
The outbreak of COVID-19 led to a record-breaking race to develop a vaccine. However, the limited vaccinecapacity creates another massive challenge: how to distribute vaccines to mitigate the near-end impact of thepandemic? In the United States in particular, the new Biden administration is launching mass vaccinationsites across the country, raising the obvious question of where to locate these clinics to maximize the effec-tiveness of the vaccination campaign. This paper tackles this question with a novel data-driven approach tooptimize COVID-19 vaccine distribution. We first augment a state-of-the-art epidemiological model, calledDELPHI, to capture the effects of vaccinations and the variability in mortality rates across age groups. Wethen integrate this predictive model into a prescriptive model to optimize the location of vaccination sitesand subsequent vaccine allocation. The model is formulated as a bilinear, non-convex optimization model.To solve it, we propose a coordinate descent algorithm that iterates between optimizing vaccine distributionand simulating the dynamics of the pandemic. As compared to benchmarks based on demographic and epi-demiological information, the proposed optimization approach increases the effectiveness of the vaccinationcampaign by an estimated 20%, saving an extra 4,000 extra lives in the United States over a three-monthperiod. The proposed solution achieves critical fairness objectives—by reducing the death toll of the pan-demic in several states without hurting others—and is highly robust to uncertainties and forecast errors—byachieving similar benefits under a vast range of perturbations.
Key words : COVID-19; Epidemiological modeling; Vaccine distribution; Non-convex optimization
1. Introduction
The outbreak of the COVID-19 pandemic has started a global race to develop vaccines, fueledby extensive investments, governmental support, and scientific breakthroughs. Thanks to theseunprecedented efforts, the scientific community delivered the good news that the whole world waseagerly awaiting. By Summer 2020, several vaccines had been developed. By the end of 2020,several vaccines got approved for emergency use and hundreds more were going under developmentand testing. Whereas vaccine development used to take years and even decades, these results rank,with no doubt, among the greatest scientific achievements (Lurie et al. 2020, Graham 2020).Unfortunately, discovering and developing a vaccine for COVID-19 was just the beginning—itwill now take months to produce, distribute, and deliver vaccines at scale. The world has quicklycome to the realization that vaccines cannot be made available immediately to everyone, andpolicy makers need to make tough decisions to pilot vaccine distribution. A global consensus hasnaturally emerged to prioritize to healthcare workers, other front line workers, and vulnerable a r X i v : . [ q - b i o . P E ] F e b ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? populations such as older people and people with comorbidities (see, e.g., National Academies ofSciences, Engineering, and Medicine 2020). Within these general principles, each jurisdiction isdesigning more detailed eligibility guidelines to distribute vaccines effectively and equitably withina population , based on demographic, clinical and geographic factors. However, a question remainsopen: how to plan vaccine distribution across populations , that is, how to allocate a limited vaccinesupply across communities, across provinces, and even across countries?In the United States, this question gained prominence in the midst of a presidential transition.In particular, the new Biden administration relies on higher extents of federal coordination invaccine distribution, as opposed to a more decentralized approach at the state level. In one of itsfirst major decisions, the administration started opening mass vaccination sites, with many moreplanned over the next few weeks. This environment raises the critical question of where to locatethese vaccination sites. Obviously, these decisions need to adhere to a number of political andfairness considerations—most notably, there must be at least one site per state. Yet, there remainsflexibility to use mass vaccination sites as a strategic lever to effectively combat the pandemic.This paper addresses this question with a novel data-driven approach, combining epidemiologicalmodeling and prescriptive analytics, to optimize the location of vaccination sites and the subsequentallocation of vaccines. To this end, we leverage a recent compartmental epidemiological modelcalled DELPHI (Differential Equations Lead to Predictions of Hospitalizations and Infections),which extends Susceptible-Exposed-Infected-Recovered (SEIR) models to capture critical drivers ofthe COVID-19 pandemic: (i) under-detection due to limited testing, (ii) governmental and societalresponse, and (iii) declining mortality rates (Li et al. 2020). The DELPHI model has been fittedfrom historical data at the country level, at the state level in the United States, and at the provincelevel in a few other countries. The DELPHI forecasts have been incorporated into the ensembleforecast from the US Center for Disease Control (2020a) and have been utilized in selecting thePhase III trial locations for the Johnson and Johnson COVID-19 vaccine. Historically, the DELPHImodel has featured excellent predictive performance, matching the number of detected cases anddeaths with high accuracy across the various waves of the pandemic.In this paper, we integrate the (predictive) DELPHI model into a (prescriptive) optimizationmodel for vaccine allocation. We first propose an extension of DELPHI, referred to as DELPHI–V,to capture the effects of vaccinations on the dynamics of the pandemic. The DELPHI–V modelalso disaggregates the dynamics of the pandemic at the subpopulation level to reflect disparitiesin mortality rates across age groups, which are critical drivers of vaccination strategies. We thenformulate an optimization model, referred to as DELPHI–V–OPT, which optimizes the vaccine ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? distribution strategy (that is, the deployment of mass vaccination sites at the strategic level, and thesubsequent allocation of vaccines at the tactical level) to minimize the death toll of the pandemic.From a technical standpoint, the DELPHI–V–OPT model relies on time discretization to embedthe system of ordinary differential equations governing the DELPHI–V dynamics into an opti-mization model. The model is formulated as a bilinear (non-convex) optimization model, due tothe SEIR dynamics at the core of DELPHI–V in which the number of new cases is driven by thenumber of susceptible and infected people. To solve it efficiently in realistic large-scale settings, wepropose a coordinate descent algorithm. Starting from a baseline solution, the algorithm iterates,until convergence, between optimizing the vaccine distribution strategy (for given dynamics of thepandemic) and simulating the dynamics of the pandemic (for a given vaccine distribution strategy).We implement the proposed model and algorithm using real-world data in the United Statesfrom the New York Times (2020), the US Census Bureau (2020), and the US Center for DiseaseControl (2021). We leverage the parameter estimates from the DELPHI model in each US state.One challenge, however, is that DELPHI estimates mortality rates in each state in each time period,while the US Center for Disease Control (2021) reports mortality rates in each age bracket. Todevelop realistic and consistent estimates for mortality rates in each state, each age group and eachtime period, we formulate another bilinear optimization model that interpolates these two piecesof information, while ensuring consistency with broader demographic information.Results suggest that the locations of vaccination sites can have a massive impact on the effec-tiveness of the vaccination campaign. As compared to several benchmarks based on demographicinformation (e.g., city and state population) and epidemiological information (e.g., case counts),our optimization approach increases the number of lives saved by the vaccines by 20%, or 4,000 livesover a three-month period in the United States. These results underscore the necessity to considerboth demographics and epidemiological dynamics when determining the locations of vaccinationsites and subsequent vaccine allocation, which is achieved by the combination of our DELPHI–Vepidemiological model and our optimization framework. In addition, the optimization approachcan ensure equity between states and across vaccination sites, thus alleviating the death toll ofthe pandemic in some states without hurting others. Finally, these benefits are highly robust tomisspecifications and fluctuations in the DELPHI parameters. Practically speaking, even thoughtactical decisions (e.g., vaccine allocation) need to be revised continuously in response to the latestinformation available throughout the vaccination campaign, strategic decisions (i.e., the locationof vaccination sites) are highly robust to noise and uncertainty.In summary, this paper makes three contributions. From a modeling standpoint, it formulates anovel optimization model for vaccine allocation, DELPHI–V–OPT, that integrates a state-of-the-art epidemiological model into an optimization model that supports vaccine distribution strategies, ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? in order to mitigate the impact of the pandemic. From a computational standpoint, it develops ascalable coordinate descent algorithm, which converges effectively and in short runtimes. From apractical standpoint, it demonstrates that optimizing the locations of mass vaccination sites cancurb the death toll of COVID-19 by a sizeable amount, thus highlighting the critical role of vaccinedistribution besides vaccine design and vaccine production in combating the pandemic. Obviously,vaccine distribution involves broad political, economic and social considerations, which lie beyondthe scope of this paper; yet, this paper can play a critical role to support ongoing mass vaccinationefforts in order to mitigate the impact of the pandemic on public health.
2. Literature review
Many pharmaceutical companies and academic institutions have explored different technologiestoward a SARS-CoV-2 vaccine (Shin et al. 2020, Florindo et al. 2020). These span (i) inactivated orlive-attenuated virus vaccines, which induce an immune response from weakened or killed pathogens(used by the Wuhan Institute of Biological Products, for instance); (ii) viral vector vaccines, whichexploit non-replicating adenoviruses to deliver an antigenic element (used by Johnson and Johnson,for instance); (iii) subunit vaccines, which use a minimal structural component of a pathogen suchas a protein (used by Clover Biopharmaceuticals, for instance); (iv) nucleic acid vaccines, whichdeliver DNA or mRNA of viral proteins (used by Pfizer and Moderna, for instance).From an operational standpoint, a vast literature studies vaccine supply chains (see Duijzeret al. 2018, Lemmens et al. 2016). A first area involves optimizing vaccine composition (Wu et al.2005, Kornish and Keeney 2008, Cho 2010, Bandi and Bertsimas 2020). A second area focuseson vaccine production to manage supply-side and demand-side uncertainty and mitigate incentivemisalignments between manufacturers and end users (Chick et al. 2008, Federgruen and Yang 2009,Arifo˘glu et al. 2012). Next, vaccine allocation optimizes the management of a vaccine stockpile(Sun et al. 2009, Mamani et al. 2013). Last, vaccine delivery optimizes inventory, distributionand dispensing operations (Jacobson et al. 2006a, Aaby et al. 2006, Dai et al. 2016). Most ofthis research focuses on predictable and repeatable epidemics, such as seasonal influenza. Forless predictable epidemics, such as pandemic influenza, advance planning interventions includestockpiling (Jacobson et al. 2006b) and anticipatory vaccination (Arinaminpathy et al. 2012).Unfortunately, these approaches are not readily applicable to a new disease such as COVID-19.Our paper deals with centralized vaccine allocation within a population. Early studies establishedthe importance of partitioning the population into risk classes (e.g., age groups) to reflect theimpact of an epidemic (Watson 1972, Elveback et al. 1976, Longini Jr et al. 1978). Emanuel andWertheimer (2006) propose a life-cycle model that prioritizes the most valuable subpopulations.Within a region, results suggest prioritizing at-risk populations (Patel et al. 2005, Chowell et al. ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites?
3. Model formulation
Our model optimizes vaccine distribution strategy. In the US context, this primarily involves thelocation of mass vaccination sites. However, optimizing these decisions requires to account forsubsequent vaccine allocation across the population, in order to further optimize and evaluate the ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites? effects of the vaccination campaign. Therefore, we refer to as vaccine distribution strategy the setof three decisions: (i) the location of mass vaccination sites, (ii) the allocation of vaccines acrossvaccination sites, and (iii) the allocation of vaccines within each sub-population.We capture the dynamics of the pandemic by means of an epidemiological model, called DELPHI,which forecasts the number of detected cases, hospitalizations and deaths in each US state (Li et al.2020). We review it briefly, and augment it to capture the effects of vaccinations—we refer to thismodel as DELPHI–V. We then embed the DELPHI–V model into a mathematical programmingmodel to optimize vaccine allocation, referred to as DELPHI–V–OPT.
DELPHI is a compartmental epidemiological model, which extends the widely used SEIR modelto account for specificities of the COVID-19 pandemic. The model is governed by a system ofordinary differential equations (ODEs) across 11 states: susceptible ( S ), exposed ( E ), infectious( I ), undetected cases who will recover ( U R ) or die ( U D ), hospitalized cases who will recover ( H R )or die ( H D ), quarantined cases who will recover ( Q R ) or die ( Q D ), recovered ( R ) and dead ( D ).DELPHI differs from most other COVID-19 forecasting models (see, e.g. Kissler et al. 2020,Perkins and Espana 2020, Rodriguez et al. 2020) by capturing three key elements of the pandemic: • Under-detection : Many cases remain undetected due to limited testing, asymptomatic car-riers, and detection errors. Ignoring them would underestimate the scale of the pandemic. TheDELPHI model captures them through the U R and U D states. • Governmental and societal response : Social distancing policies limit the spread of thevirus. Ignoring them would overestimate the scale of the pandemic. However, if restrictions are liftedprematurely, a resurgence may occur. We define a governmental and societal response function γ ( t ), which modulates the infection rate and is parameterized as follows: γ ( t ) = 1 + 2 π arctan (cid:18) − ( t − t int ) ω (cid:19) + c exp (cid:18) − ( t − t jump ) σ (cid:19) . (1)This parameterization defines four phases (Figure 1). In Phase I, most activities continue normally.In Phase II, the infection rate declines sharply as policies get implemented. The parameters t int and ω can be interpreted as the start time and strength of this response. In Phase III, the declinereaches saturation. The epidemic then experiences a resurgence of magnitude c in Phase IV, dueto relaxations in governmental and social restrictions. This is counteracted at time t jump , whenrestrictions are re-implemented, with σ controlling the duration of this second wave. DELPHI is also applied to each country and to other provinces, but this paper focuses on US states. ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites? Figure 1 Governmental and societal response function γ ( t ) ( ω = 5 , t int = 10 , c = 1 , t jump = 25 and σ = 2 ). • Declining mortality rates : The mortality rate of COVID-19 has been declining throughthe pandemic, due to a better detection of mild cases, enhanced care for COVID-19 patients, andother factors. We model the mortality rate as a monotonically decreasing function of time: m ( t ) = ( m − m min ) (cid:18) π arctan ( − r m t ) (cid:19) + m min , (2)where m is the initial mortality rate, m min is the minimum mortality rate and r m is a decay rate.Ultimately, DELPHI involves 16 parameters that define the transition rates between the 11states. We calibrate seven of them from a database on clinical outcomes (Bertsimas et al. 2020).Using non-linear optimization, we estimate the other 9 parameters from historical data on thenumber of cases and deaths in each region. We refer to Li et al. (2020) for details.Since its inception in March 2020, DELPHI has been extensively tested and validated againstreal-world data. Figure 2 reports the historical performance of the model in the United States,during the first wave in the Spring of 2020 and the second wave in the Fall of 2020. As the figureshows, the model has been predicting the magnitude of the pandemic with high accuracy up toone month in advance; for instance, as early as April 3, 2020, the model was predicting 1.2–1.4million cases in the United States by early May, a prediction that became quite accurate a monthlater (Figure 2a). Obviously, subsequent forecasts, by leveraging more up-to-date information, wereable to refine these estimates. As a result, the DELPHI model was incorporated into the ensembleforecast from the US Center for Disease Control (2020a). During the second wave of the pandemic,DELPHI continued to exhibit strong predictive performance, with a mean average percentage erroramong the lowest of the CDC ensemble forecast (Figure 2b). ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? (a) US First-wave predictions (b) Second-wave performance Figure 2 Historical performance of the DELPHI predictions in the United States.
We now augment the DELPHI model to capture two key aspects of vaccinations:1.
Disparate impacts of the disease across risk classes . Age is one of the primary drivers ofmortality (Guan et al. 2020, Goyal et al. 2020, Petrilli et al. 2020). The US Center for DiseaseControl (2020b) reports that the mortality rate among Americans aged 70 and over is twoorders of magnitude greater than for those aged 30 and under. We partition the populationinto risk classes , defined as homogeneous groups with comparable health characteristics. Weconsider age-based risk classes in our experiments, but other categorizations could be used(e.g., based on comorbidities). Accordingly, we replicate the 11 model states for each risk class.2.
Impact of vaccinations on the dynamics of the pandemic . A fraction of vaccinated peoplewill be immune to the disease (based on the vaccine’s effectiveness). Clinical trials suggestthat early-approved vaccines prevent mortality but not necessarily infections. Therefore, weassume conservatively that all vaccinated people can still transmit the disease. We relax thisassumption later on, to show the robustness of our results when a fraction of vaccinatedpeople become fully immune to the disease. We create four new model states: susceptible andvaccinated ( S (cid:48) ), exposed and vaccinated ( E (cid:48) ), infected and vaccinated ( I (cid:48) ), and immune ( M ).Figure 3 shows a simplified flow diagram of the DELPHI–V model, with two risk classes (indexedby k = 1 , ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? optimization model. Accordingly, we denote the states of undetected, hospitalized and quarantinedpeople who will die from the disease by U , H and Q (as opposed to U D , H D and Q D ). VaccinatedAge group 1Age group 2 S E I U H Q D S E I U H Q D S (cid:48) E (cid:48) I (cid:48) MβV βV Figure 3 Simplified flow diagram of the DELPHI–V model.
Let V k ( t ) denote the population mass from risk class k ∈ K that gets vaccinated at time t ∈ T , andlet β ∈ (0 ,
1] denote the vaccine’s effectiveness. For simplicity, we make three assumptions. First,the effects of vaccines are instantaneous (relaxing this assumption, although straightforward, wouldmerely induce a time lag into the system, without significantly impacting the vaccine distributionstrategy). Second, the vaccine has no effect when it fails to immunize the patient (i.e., no partialbenefit and no side effect). Third, we consider single-dose vaccines (double-dose vaccines wouldagain induce a time lag without impacting the vaccine distribution strategy). Then, at each time t ∈ T , a mass βV k ( t ) of people transitions from the susceptible state S k to the state S (cid:48) , and theremaining mass (1 − β ) V k ( t ) remains in the susceptible state. People in the S (cid:48) state can becomeexposed and infected, but then become immune to the disease (as opposed to having a positiveprobability of dying from it). Note that infections are driven by the total mass of infected people,across all risk classes and vaccinated people. All other transitions shown in Figure 3 are consistentwith the DELPHI model.The DELPHI–V model is governed by the following ODE system:d S k d t = − βV k ( t ) − αγ ( t ) ( S k ( t ) − βV k ( t )) (cid:32) K (cid:88) l =1 I l ( t ) + I (cid:48) ( t ) (cid:33) (3)d S (cid:48) d t = + β K (cid:88) κ =1 V κ ( t ) − αγ ( t ) (cid:32) S (cid:48) ( t ) + β K (cid:88) κ =1 V κ ( t ) (cid:33) (cid:32) K (cid:88) l =1 I l ( t ) + I (cid:48) ( t ) (cid:33) (4)d E k d t = αγ ( t ) ( S k ( t ) − βV k ( t )) (cid:32) K (cid:88) l =1 I l ( t ) + I V ( t ) (cid:33) − r I E k ( t ) (5) ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? d E (cid:48) d t = αγ ( t ) (cid:32) S (cid:48) ( t ) + β K (cid:88) κ =1 V κ ( t ) (cid:33) (cid:32) K (cid:88) l =1 I l ( t ) + I (cid:48) ( t ) (cid:33) − r I E (cid:48) ( t ) (6)d I k d t = r I E k ( t ) − r d I k ( t ) (7)d I (cid:48) d t = r I E (cid:48) ( t ) − r d I (cid:48) ( t ) (8)d U k d t = r Uk ( t ) I k ( t ) − r D U k ( t ) (9)d H k d t = r Hk ( t ) I k ( t ) − r D H k ( t ) (10)d Q k d t = r Qk ( t ) I k ( t ) − r D Q k ( t ) (11)d D k d t = r D ( U k ( t ) + H k ( t ) + Q k ( t )) (12)d M d t = r d I (cid:48) ( t ) , (13)where:– α is the nominal infection rate;– γ ( t ) : R + → R + is the governmental and societal response function (Figure 1);– r I , r d , r D , are the progression rate, the detection rate, and the death rate;– r Uk ( t ), r Hk ( t ), and r Qk ( t ) capture the detection, hospitalization and death rates, accounting forthe probabilities of detection and hospitalization and the mortality rate (Equation (2)). Theirdependency on t and k reflect disparities over time and across risk classes.As noted earlier, the dynamics of exposure and infection depend on the total number of infectedpeople (across risk classes and vaccinated/non-vaccinated people), as opposed to the number ofinfected people in a given risk class. DELPHI–V captures these interdependencies—indicated bythe red rectangle in Figure 3 and the terms (cid:80) Kl =1 I l ( t ) + I (cid:48) ( t ) in Equations (3)–(5).Given initial conditions, the ODE equations uniquely determine the evolution of this systemover time—for a given vaccine allocation reflected in the variable V . Next, we optimize the vaccinedistribution strategy to minimize the overall impact of the pandemic—estimated by DELPHI–V. The DELPHI–V–OPT model takes as inputs epidemiological information (estimated from theDELPHI–V model), information on the vaccine (including vaccine effectiveness and vaccine bud-get), and demographic information in the United States (e.g., major cities, distance across counties,population per county). It optimizes the vaccine distribution strategy, including the location ofmass vaccination sites and the subsequent allocation of vaccines. It is formulated as a tri-objectivemodel, to minimize (i) the death toll of the pandemic, (ii) the number of exposed people in thetermination period, and (iii) the distance between vaccination sites and population centers. The ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites? main public health objective is obviously death minimization, so the first objective componentis heavily prioritized. However, just considering the number of deaths could result in a waste ofvaccines near the end of the planning horizon, as individuals infected in the final periods would nothave time to flow to the death state in the epidemiological model. Therefore, the second componentof the objective minimizes the number of infections. The last component minimizes geographicdisparities. In addition, the model incorporates other equity consideration by means of fairnessconstraints.We proceed by time discretization to formulate the optimization model and retain tractability.This reduces to solving the system of ODE equations given in Equations (3)–(13) by a forwarddifference scheme. We denote by ∆ t the discretization unit (e.g., 1 day).Formally, we define the following sets, input parameters, and decision variables. Sets L = set of population centers in the United States (e.g., counties), { , · · · , L }M = set of regions in the United States (e.g., 50 states plus the District of Columbia), { , · · · , m }N = set of candidate vaccination sites, { , · · · , n }N j = subset of candidate vaccination sites located in state j ∈ MN l = subset of candidate vaccination sites located in the same state as county l ∈ LK = set of risk classes in the population (e.g., age groups), { , · · · , K }T = set of discretized time periods (e.g., days), { , · · · , T } Parameters
P op l = total population in center l ∈ L ∆ li = distance from population center l ∈ L to candidate vaccination site i ∈ N N = number of vaccination sites to be deployed across the country B t = number of vaccines available at time t ∈ T β = effectiveness of vaccines r I = progression rate of the disease (DELPHI parameter) r d = detection rate of the disease (DELPHI parameter) r D = death rate of the disease (DELPHI parameter) r Ujkt = transition rate from I to U , capturing mortality ratein region j ∈ M at time t ∈ T for risk class k ∈ K (DELPHI parameter) r Hjkt = transition rate from I to H , capturing mortality ratein region j ∈ M at time t ∈ T for risk class k ∈ K (DELPHI parameter) r Qjkt = transition rate from I to Q , capturing mortality ratein region j ∈ M at time t ∈ T for risk class k ∈ K (DELPHI parameter) r Hjkt = transition rate from I to H , capturing mortality rate α j = nominal infection rate in region j ∈ M (DELPHI parameter) γ jt = governmental and societal response in region j ∈ M at time t ∈ T (DELPHI parameter) ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? Note that the parameters r U , r H and r Q are defined for each region, risk class and time pe-riod, reflecting underlying variations in mortality rates. In contrast, the parameters r I , r d and r D are treated as uniform characteristics of the disease. In reality, these parameters may varyacross risk classes; for instance, the serological estimates from the US Center for Disease Control(2020a) suggest different prevalence of the disease across age groups. We test this hypothesis inour experiments, to verify the robustness of our results to the uniform infection rate assumption.We also assume a single vaccine effectiveness value β . In theory, vaccine effectiveness mightalso vary across risk classes. More importantly, there are now several vaccines available, each withdifferent clinical characteristics. Ideally, we could introduce an additional set to capture vaccineheterogeneity, and let the vaccine effectiveness vary across vaccine types. This approach however,may be somewhat impractical in practice, as it may be difficult to strategically allocate differentvaccines to different populations based on vaccine effectiveness. For equity, we therefore assumeconservatively that the mix of vaccines remains identical across vaccination sites. Under this re-striction, the mix of vaccines can be reduced to a representative vaccine with average effectiveness. Primary decision variables x i = (cid:40) i ∈ N is selected0 otherwise C it : number of vaccines distributed to site i ∈ N at time t ∈ T W li = (cid:40) l ∈ L is assigned to vaccination site i ∈ N S jkt : number of eligible people to region j ∈ M in risk class k ∈ K at time t ∈ T V jkt : number of vaccines allocated to region j ∈ M in risk class k ∈ K at time t ∈ T To track the impact of vaccine allocation on the resulting dynamics of the pandemic, we createindirect variables, corresponding to all the states in the DELPHI–V model shown in Figure 3.The vaccine distribution problem, referred to as ( P ), is then formulated in Equations (14)–(39). min M (cid:88) j =1 K (cid:88) k =1 ( D jkT + H jkT + Q jkT ) + λ E M (cid:88) j =1 K (cid:88) k =1 E jkT + λ D L (cid:88) l =1 n (cid:88) i =1 P op l ∆ li W li (14)s.t. n (cid:88) i =1 x i = N (15) (cid:88) i ∈N j x i ≥ , ∀ j ∈ M (16) W li ≤ x i , ∀ l ∈ L , i ∈ N (17) (cid:88) i ∈N l W li = 1 , ∀ l ∈ L (18) ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? n (cid:88) i =1 C it ≤ B t , ∀ t ∈ T (19) C it ≤ B t x i , ∀ i ∈ N , t ∈ T (20) K (cid:88) k =1 V jkt ≤ (cid:88) i ∈N j C it , ∀ j ∈ M , t ∈ T (21)¯ S j,k,t +1 ≤ ¯ S jkt − βV jkt − ( S jkt − S j,k,t +1 ) , ∀ j ∈ M , k ∈ K , t ∈ T (22) V jkt ≤ ¯ S jkt , ∀ j ∈ M , k ∈ K , t ∈ T (23) | C i,t +1 − C it | ≤ θ S C it , ∀ i ∈ N , t ∈ T (24) (cid:80) l ∈L j P op l (cid:80) mj (cid:48) =1 (cid:80) l ∈L j (cid:48) P op l N − Θ L ≤ (cid:88) i ∈N j x i ≤ (cid:80) l ∈L j P op l (cid:80) mj (cid:48) =1 (cid:80) l ∈L j (cid:48) P op l N + Θ L , ∀ j ∈ M (25) B t N (1 + θ V ) x i ≤ C it ≤ B t N (1 + θ V ) x i , ∀ i ∈ N , t ∈ T (26) (cid:88) i ∈N j C it ≤ (cid:32) (cid:80) l ∈L j P op l (cid:80) mj (cid:48) =1 (cid:80) l ∈L j (cid:48) P op l + θ P (cid:33) B t , ∀ j ∈ M , t ∈ T (27) S j,k,t +1 ≥ S jkt − βV jkt − α j γjt ( S jkt − βV jkt ) (cid:32) K (cid:88) κ =1 I jκt + I (cid:48) jt (cid:33) ∆ t, ∀ j ∈ M , k ∈ K , t ∈ T (28) S (cid:48) j,t +1 ≥ S (cid:48) jt + β K (cid:88) κ =1 V jκt − α j γ jt (cid:32) S (cid:48) jt + β K (cid:88) κ =1 V jκt (cid:33) (cid:32) K (cid:88) κ =1 I jκt + I (cid:48) jt (cid:33) ∆ t, ∀ j ∈ M , t ∈ T (29) E j,k,t +1 ≥ E jkt + (cid:32) α j γ jt ( S jkt − βV jkt ) (cid:32) K (cid:88) κ =1 I jκt + I (cid:48) jt (cid:33) − r I E jkt (cid:33) ∆ t, ∀ j ∈ M , k ∈ K , t ∈ T (30) E (cid:48) j,t +1 ≥ E (cid:48) jt + (cid:32) α j γ jt (cid:32) S (cid:48) jt + β K (cid:88) κ =1 V jκt (cid:33) (cid:32) K (cid:88) κ =1 I jκt + I (cid:48) jt (cid:33) − r I E (cid:48) jt (cid:33) ∆ t, ∀ j ∈ M , t ∈ T (31) I j,k,t +1 ≥ I jkt + (cid:0) r I E jkt − r d I jkt (cid:1) ∆ t, ∀ j ∈ M , k ∈ K , t ∈ T (32) I (cid:48) j,t +1 ≥ I (cid:48) jt + (cid:0) r I E (cid:48) jt − r d I (cid:48) jt (cid:1) ∆ t, ∀ j ∈ M , t ∈ T (33) U j,k,t +1 ≥ U jkt + (cid:0) r Ujkt I jkt − r D U jkt (cid:1) ∆ t, ∀ j ∈ M , k ∈ K , t ∈ T (34) H j,k,t +1 ≥ H jkt + (cid:0) r Hjkt I jkt − r D H jkt (cid:1) ∆ t, ∀ j ∈ M , k ∈ K , t ∈ T (35) Q j,k,t +1 ≥ Q jkt + (cid:0) r Qjkt I jkt − r D Q jkt (cid:1) ∆ t, ∀ j ∈ M , k ∈ K , t ∈ T (36) D j,k,t +1 = D jkt + r D ( U jkt + H jkt + Q jkt ) ∆ t, ∀ j ∈ M , k ∈ K , t ∈ T (37) M j,t +1 = M jkt + r d I (cid:48) jt , ∀ j ∈ M , t ∈ T (38) x , W binary , C , V , S , S (cid:48) , ¯S , E , E (cid:48) , I , I (cid:48) , U , H , Q , D , M ≥ . (39) Equation (14) formalizes the three objectives of the model. The first term corresponds to ourprimary objective of minimizing the number of deaths over the planning horizon, across all regionsand risk classes. This number includes people in the absorbing state D , as well as the transient states H and Q (we ignore undetected deaths). The next terms minimize, as lower-priority objectives, thenumber of exposed people at the end of the horizon and the distance to the vaccination sites. Thehyperparameters λ E and λ D are set to small values to prioritize the death-minimization objective.Next, the constraints capture practical considerations surrounding vaccine distribution: ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? – Number of vaccination sites : We impose a total budget of N vaccination sites (Equa-tion (15)). For obvious reasons, there needs to be at least one site in every state (Equation (16)).We consider N = 100 in our experiments, which leaves flexibility to strategically deploy 49 sites.– Assignment : Equation (17) ensures that people get assigned to vaccination sites that havebeen selected. Equation (18) assigns each population center to exactly one site, in the same state.These assignment constraints are used to compute the distance term in the objective function.–
Inter-regional vaccine capacity : Due to restrictions in vaccine manufacturing and distri-bution networks, a limited number of vaccines can be allocated in each time period. Equation (19)ensures that the total number of vaccines allocated lies within the available budget in each period.–
Consistency : Equation (20) ensures that vaccines only get distributed to selected sites. Simi-larly, Equation (21) ensures that the number of people vaccinated in each state (across risk classes)does not exceed the number of vaccines allocated that state. This constraint involves two assump-tions. A first, conservative assumption is that people can only get vaccinated in the state that theylive in, which is required in practice for traceability purposes. Another, optimistic assumption isthat the vaccine allocation constraint applies to each state, as opposed to each vaccination site.In other words, the model assumes vaccines can be reallocated between vaccination sites within astate, thus maintaining a degree of freedom in intra-state vaccine distribution.–
Eligibility : We prevent people from being vaccinated twice: a patient who has been vaccinatedbut remains susceptible cannot be vaccinated again. Equation (22) defines the number eligiblepeople as the previous number of eligible people minus the number of people for whom the vaccinewas effective and the number of people who got exposed to the disease. Equation (22) then ensuresthat the number of vaccinated people lies below the number of eligible people.–
Smoothness : Large fluctuations in the number of vaccines allocated to each region from dayto day would likely cause problems from a supply chain management perspective—both to deliverand to administer the vaccines. Equation (24) ensures that such fluctuations remain minimal. Thehyperparameter θ S controls the trade-off between efficiency and smoothness.– Fairness : To be politically and socially viable, vaccine distribution must not neglect anyregion, even if it is not a virus “hot spot”. This also enhances the robustness of the solution, giventhat inter-regional transmission can occur in practice. Equation (25) promotes inter-state fairnessat the strategic level, by ensuring that the fraction of vaccination sites in each state does notdeviate too much from its population share. Equation (26) promotes inter-site fairness, by ensuringthat vaccine distribution across sites does not deviate too much from uniform distribution. Finally,Equation (27) promotes inter-state fairness at the tactical level, by ensuring that no state receives afraction of vaccines that exceeds its population share by a wide margin. The hyperparameters Θ L , θ V and θ P control the trade-off between efficiency and fairness. As the results will show, even tightfairness constraints leave critical flexibility when locating vaccination sites and allocating vaccines. ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? – DELPHI–V dynamics : Equations (28)–(38) capture the dynamics of the DELPHI–V modelin a discretized time space (Equations (3)–(13)).–
Domain of definition : Equation (39) defines the domain of each variable.
Model Structure
Problem ( P ) is a non-linear program, due to the bilinear terms in Equations (28)–(31), whichreflect the fact that the number of new infections result from the interactions between susceptibleand infected populations—a key characteristic of all SEIR-based compartmental models. Thesebilinear terms result in non-convex constraints, thus in a highly challenging optimization model.The latest Gurobi 9.0 release includes a solver for non-convex quadratic problems (Gurobi Op-timization 2020). Yet, general-purpose technologies are limited to small-scale instances. In oursetting, Problem ( P ) includes 2 m ( K + 1) T non-convex constraints each involving 2( K + 1) bilinearterms, for a total of 4 mT ( K + 1) bilinear terms. A realistically-sized problem with M = 51 (50 USstates plus Washington, D.C.), K = 6 (6 age groups) and T = 90 (a 3 month planning horizon withdaily discretization) would result in nearly 900,000 bilinear terms. Problem ( P ) remains intractablewith existing commercial solvers, motivating the development of a tailored algorithm.
4. Solution algorithm
We propose an iterative coordinate descent algorithm to solve Problem ( P ) in short computationaltimes—consistent with practical requirements. We describe the algorithm in this section. We alsopresent three baselines replicating reasonable strategies that could be implemented in the absenceof our data-driven optimization model. These baselines are used for two purposes: (i) to providean initial feasible solution in the coordinate descent algorithm, and (ii) as benchmarks to evaluatethe benefits of the data-driven optimization approach proposed in this paper. Our algorithm relies on two key observations: 1) aside from Equations (28)–(31), the objectivefunction and all other constraints in ( P ) are linear, and 2) given a fixed vaccine distribution strategy,the discretized DELPHI–V model can be solved efficiently. Therefore, we proceed by coordinatedescent, alternating between two modules: one that optimizes the vaccination distribution strategygiven the infection dynamics, and one that simulates the bilinear dynamics of the pandemic for agiven vaccination distribution strategy. The optimization part reduces to a linear program, whichcan be solved very efficiently. Using the resultant vaccine distribution, the simulation part re-estimates the infected population under bilinear dynamics, using a forward discretization scheme.Specifically, the two modules are defined as follows: ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? Simulate : Based on a vaccine allocation solution V , we compute the DELPHI–V dynamicsfrom t = 0 to t = T (Section 3.2) by solving the ODE system (Equations (3)–(13)) using aforward difference scheme in a discretized time space. This terminates in O ( M KT ) operations.We denote the total infected population (across all risk classes and vaccinated people) in region j at time t ∈ T by (cid:98) I jt = (cid:80) Kk =1 I jkt + I (cid:48) jt . We refer to this procedure as Simulate ( x , C , W , V ).2. Optimize : Given the infectious population estimates (cid:98) I , we can approximate Equation (28)–(31) by the following linear constraints. The problem can then be efficiently solved as a linearprogramming model. We refer to this module as Optimize (cid:16)(cid:98) I (cid:17) . S j,k,t +1 ≥ S jkt − βV jkt − α j γ jt ( S jkt − βV jkt ) (cid:98) I jt ∆ t, ∀ j ∈ M , k ∈ K , t ∈ T S (cid:48) j,t +1 ≥ S (cid:48) jt + β K (cid:88) κ =1 V jκt − α j γ jt (cid:32) S (cid:48) jt + β K (cid:88) κ =1 V jκt (cid:33) (cid:98) I jt ∆ t, ∀ j ∈ M , t ∈ T E j,k,t +1 ≥ E jkt + (cid:16) α j γ jt ( S jkt − βV jkt ) (cid:98) I jt − r I E jkt (cid:17) ∆ t, ∀ j ∈ M , k ∈ K , t ∈ T E (cid:48) j,t +1 ≥ E (cid:48) jt + (cid:32) α j γ jt (cid:32) S (cid:48) jt + β K (cid:88) κ =1 V jκt (cid:33) (cid:98) I jt − r I E (cid:48) jt (cid:33) ∆ t, ∀ j ∈ M , t ∈ T We iterate between the
Simulate and
Optimize modules, until convergence. Specifically, thealgorithm terminates when the variation in the objective function value remains minimal from oneiteration to the next. The pseudocode summarizing this approach is presented in Algorithm 1. Weturn next to the generation of an initial feasible solution.
We propose three simple and interpretable baselines for generating a feasible solution to ( P ). Bydesign, these baselines are heuristics that solely rely on the inputs of the optimization models, asopposed to requiring the full model and algorithm developed in this paper. Top-cities baseline : This approach prioritizes cities based on population. Specifically, it deploysvaccination sites in the most populous cities, while accounting for the constraint that each statemust have at least one center (Equation (16)). Subsequently, it allocates an equal fraction of thedaily vaccine budget to each vaccination site: we fix the variables x and C , run the model tooptimize vaccine allocation within each state, and estimate the resulting number of deaths. Thisbaseline corresponds to a city-level approach based on demographic information alone. Population-based baseline : Under this approach, the number of vaccination sites deployed ineach state is based on the state’s population share. This is formulated as follows, where y j is adecision variable denoting the number of vaccination sites in state j .min m (cid:88) j =1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:80) l ∈L j P op l (cid:80) mj (cid:48) =1 (cid:80) l ∈L j (cid:48) P op l N − y j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , s.t. m (cid:88) j =1 y j = 100, y j ≥ , ∀ j ∈ M , y integer. ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? Algorithm 1:
Coordinate descent algorithm for DELPHI–V–OPT (Problem ( P )). input : Prescriptive DELPHI–V–OPT data; termination tolerance ε Initialization: i ← x ( i ) ← , , C ( i ) ← , , W ( i ) ← , , V i ← , I i ← , Z i ← ∞ i ← i + 1 (cid:0) x ( i ) , C ( i ) , W ( i ) , V ( i ) , I ( i ) , Z ( i ) (cid:1) ← GenerateFeasibleSolution while (cid:12)(cid:12) Z ( i ) − Z ( i − (cid:12)(cid:12) /Z ( i − > ε (cid:48) do i ← i + 1Run Simulate (cid:0) x ( i − , C ( i − , W ( i − , V ( i − (cid:1) . Update (cid:98) I ( i ) ← (cid:98) I , where (cid:98) I is the output of Simulate (cid:0) x ( i − , C ( i − , W ( i − , V ( i − (cid:1) .Run Optimize (cid:16)(cid:98) I ( i ) (cid:17) . Update (cid:0) x ( i ) , C ( i ) , W ( i ) , V ( i ) (cid:1) ← ( x , C , W , V ), where ( x , C , W , V )is the output of Optimize (cid:16)(cid:98) I ( i ) (cid:17) . Update the objective function: Z ( i ) ← M (cid:88) j =1 K (cid:88) k =1 ( D jkT + H jkT + Q jkT ) + λ E M (cid:88) j =1 K (cid:88) k =1 E jkT + λ D L (cid:88) l =1 n (cid:88) i =1 P op l ∆ li W li endoutput: Vaccine distribution strategy (cid:0) x ( i ) , C ( i ) , W ( i ) , V ( i ) (cid:1) We then solve DELPHI–V–OPT, while fixing the aggregate number of vaccination sites per state,i.e., (cid:80) i ∈N j x i = y j , and assuming equal allocation of vaccines across vaccination sites, i.e., C it = N B t , ∀ t ∈ T . The model allocates vaccines within each state and estimates the number of deaths.This baseline corresponds to a state-level approach based on demographic information alone. Case-based baseline : Under this approach, the number of vaccination sites deployed in eachstate is based on the number of COVID-19 cases at the beginning of the planning horizon. This isformulated as follows, where y j is a decision variable denoting the number of vaccination sites instate j and Cases j denotes the case count in state j .min m (cid:88) j =1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Cases j (cid:80) mj (cid:48) =1 Cases j (cid:48) N − y j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , s.t. m (cid:88) j =1 y j = 100, y j ≥ , ∀ j ∈ M , y integer.We then proceed as with the population-based baseline, by fixing the number of vaccination sitesper state, assuming equal vaccine allocation across sites, and re-solving the model. This baselinecorresponds to a state-level approach based on epidemiological information alone.By design, these baselines satisfy all constraints of Problem ( P ), and thus provide valid initializa-tions into our coordinate descent algorithm. They also provide sensible and equitable benchmarksbased on readily-available demographic information (e.g., census data) and epidemiological in-formation (e.g., case counts), hence easily implementable. Comparisons between our optimizedsolution and these benchmarks thus estimate the benefits of vaccine distribution optimization. ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites?
5. Experimental setup
We implement the proposed model and algorithm in the United States. We select N = 100 vacci-nation sites out of the 500 most populous cities in the United States as candidate locations (set N ) We define the set M as 51 “states” (the 50 states plus the District of Columbia) and the set L as the 3,006 counties. We define six risk classes based on six relatively coarse age groups: 0-9years, 10-49 years, 50-59 years, 60-69 years, 70-79 years, and 80 years and above. These simplifiedrisk classes facilitate the practical implementation of the solution while capturing broad trends inmortality rates. We define the time horizon T as the three-month period from February to April2021, consistently with the ongoing planning horizon of the US federal government. We calibrate the model using multiple data sources (Figure 4). First, we estimate the parametersof the DELPHI model (without vaccinations) independently for each state, using historical data oncases and deaths from the New York Times (2020). We obtain a granular population breakdownby age for each state from the US Census Bureau (2020). We then run DELPHI (still, withoutvaccinations) to derive the initial susceptible, exposed and infected populations (on January 30,2021), which we distribute among the risk classes proportionally to their size.The next input is an estimate of the mortality rate per region, risk class and time period. We makeuse of the data from the US Center for Disease Control (2021), which report the total number ofconfirmed COVID-19 cases, hospitalizations and deaths by age group in the United States until theend of January 2021. In contrast, the DELPHI model fits a time-dependent mortality rate withineach region. To the best of our knowledge, there exists no other data source for nationwide casesand deaths by age group. This leaves a discrepancy between the time-independent estimates at therisk class level from the CDC, and the time-varying estimates at the region-level from DELPHI.To reconcile these data, we employ an optimization procedure that interpolates the mortality rateper region, risk class and time period. We present this approach in the next section.
Our procedure to estimate mortality rates starts from two sets of inputs: • DELPHI predictions:
Let (cid:98) X jt denote the estimated number of new detected cases in region j ∈ M and time period t ∈ T . Let (cid:98) D dj,t + τ j denote the number of deaths, where τ j is the mediandeath time after detection in region j . These quantities are aggregated across risk classes. • CDC data:
Let X CDC k and D CDC k denote the number of cases and deaths for risk class k ∈ K .These quantities are aggregated across regions and time periods.We define the reference mortality rate m jk of risk class k ∈ K based in region j as follows: m jk = D CDC k X CDC k (cid:32) (cid:80) Tt =1 (cid:98) D dj,t + τ j (cid:80) Tt =1 (cid:98) X t (cid:33) (cid:32) (cid:80) Kl =1 X CDC l (cid:80) Kl =1 D CDC l (cid:33) (40) ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? NYT data US census data CDC dataDELPHI parameters Population by re-gion and risk classMortality rateestimationTrained DELPHI(without vaccinations)Initial conditions(from DELPHI) Mortality rates byregion and risk classPrescriptive model:DELPHI–V–OPTVaccine distri-bution strategy
Figure 4 Simulation environment: raw data (blue), processed data (red), models (black) and outputs (green).
By design, this expression preserves the ratio of mortality rates between different risk classes fromthe CDC data, while correcting the mean reference mortality rate in each region across the planninghorizon to be in line with the DELPHI projections.We then estimate the mortality rate for each region j ∈ M , risk class k ∈ K , and time period t ∈ T , denoted by m jkt . We also introduce additional decision variables X jkt and D j,k,t + τ j , reflectingthe number of detected cases and the number of future deaths in region j assigned to risk class k ∈ K and time period t ∈ T . Given that the fitting procedure is done separately in each region j ∈ M , we decouple the problem at the region level—thereby considerably reducing the size of eachproblem instance. Specifically, we formulate the optimization problem given in Equations (41)–(48).for all j ∈ M : min T (cid:88) t =1 K (cid:88) k =1 (cid:34)(cid:18) m jkt − m jk m jk (cid:19) + λ M (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X jkt − p jk (cid:98) X jt p jk (cid:98) X jt (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:35) (41) ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? s.t. K (cid:88) k =1 D j,k,t + τ j = (cid:98) D j,t + τ j , ∀ t ∈ T (42) K (cid:88) k =1 X jkt = (cid:98) X jt , ∀ t ∈ T (43) m jkt · X jkt = D j,k,t + τ j , ∀ k ∈ K , ∀ t ∈ T , (44) m jkt ≥ m j,k,t +1 , ∀ k ∈ K , ∀ t ∈ T (45) m jkt ≥ m jlt , ∀ k, l ∈ K : m k ≥ m l , (46) m jkt ≤ , ∀ k ∈ K , ∀ t ∈ T (47) m , C , D ≥ . (48)The first term in Equation (41) minimizes the squared relative error between the mortality rateestimates and their reference values (Equation (40)). The second term is a regularization penaltythat minimizes deviations between the proportion of detected cases and the proportion p jk of thepopulation, in each risk class. The parameter λ M trades off these two objectives (we use λ M =0 . Table 1 Calibrated monthly mortality rates, averaged by risk class and over all states.
Our DELPHI–V–OPT model relies on a forward difference scheme to simulate the dynamics ofthe pandemic, for any vaccine allocation. If the discretization is too coarse, the algorithm willintroduce truncation errors. If, however, the discretization is too granular, computational times ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites? will be prohibitively long. After extensive experimentations, we found ∆ t = 1 day to yield the bestcompromise. It is also practical choice as it yields a day-by-day plan.The coordinate descent scheme terminates when the change in the objective function value liesbelow a pre-determined threshold. We set the termination tolerance to ε = 0 . λ E , λ V , θ S , Θ L , θ V and θ P to balanceefficiency with practical and fairness considerations, and perform sensitivity analyses on theseparameters. We consider a baseline vaccine effectiveness of 90% (in line with the values from thefirst vaccine approvals) and a baseline budget of 1 million vaccines per day (which can be viewedas 1 million single-dose vaccines, 2 million double-dose vaccines or a combination thereof). Giventhe strong underlying uncertainty, we perform sensitivity analyses to demonstrate the robustnessof the benefits of our optimization outputs to vaccine effectiveness and vaccine budget.All optimization models are implemented in Gurobi 9.0, with 2.3GHz processor and 4 cores. Weuse a barrier method to solve each linear program, with a barrier convergence tolerance of 10 − .
6. Experimental results
We now present results obtained with the modeling and algorithmic framework developed in thispaper. Our main focus is on the location of the 100 vaccination sites—decisions that have to bemade immediately—as opposed to vaccine allocation—decisions that can be revisited continuouslyover time as more information becomes available. We first evaluate the benefits of the optimizedsolution against baseline approaches, and then establish the robustness of these benefits.
Figure 5 shows the heatmap of vaccination sites across all 51 states (the exact list of proposedlocations is reported in the appendix). Recall that, by design, the optimized solution deploys atleast one vaccination site per state (Equation (16)). In addition, one of the fairness constraints(Equation (25) anchors the number of vaccination sites per state to its population share. Otherthan that, the vaccination sites are deployed strategically to curb the spread of the pandemic.As a result, the optimized number of vaccination sites varies significantly per state, as a functionof underlying demographic and epidemiological factors. For instance, the four largest states bypopulation (California, Texas, Florida and New York) receive the most vaccination sites (6 to 10each). In contrast, the smaller states receive only one vaccination site. Ultimately, the heatmapsuggests that the optimized solution targets large population centers and some of the hot spots ofthe pandemic, while ensuring equity nationwide.Then, could we have achieved a similar solution with some of our benchmarks (Section 4.2), whichdetermine the locations of the vaccination sites based on demographic and epidemiological databut do not make use of our epidemiological and optimization models? To investigate this question, ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites?
22 18 1 111621 2 221 12 11 121 121 12 121 621 322 2 111 3101 132 121
Figure 5 Number of centers per state in the proposed solution. we evaluate the dynamics of the pandemic under each of the three baselines. When it comes to theoptimization, we derive three solutions: (i) an “optimized locations” solution, which optimizes thesite locations decisions freely (without fairness constraints) but then allocates vaccines uniformlyacross the 100 selected vaccination sites; (ii) a “fully optimized” solution, which optimizes the sitelocations decisions and vaccine allocation decisions freely (without fairness constraints); and (iii) a“proposed” solution, which optimizes the site locations decisions and vaccine allocation decisionswith fairness constraints. For each solution, Table 2 reports the number of lives saved by thevaccination campaign, as compared to a solution without vaccinations.Let us start with the main observation: the optimized solutions provides significant benefits,as compared to all benchmarks. Comparing the “optimized locations” solution to the top-citiesbenchmark, we find that optimizing locations alone can save around 4,500 lives over the three-month period under consideration, enhancing the impact of the vaccination campaign by 24%.Moving to the “fully optimized” solution, we note that jointly optimizing locations and vaccineallocation achieves further benefits, by saving an extra 2,000 lives and increasing the benefitsover the benchmark to 35%. These results underscore the potential of the proposed optimizationapproach, which leverages available vaccines strategically to target the regions that need themmost. Another observation is that determining the locations of vaccination sites does not achieveall potential benefits of the vaccination campaign, underscoring the role of downstream vaccineallocation as an extra lever to combat the pandemic. ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites? Table 2 Comparison of optimized solutions and benchmarks.
Method Solution Site locations Vaccine distribution Saved lives19,045Top-cities Most populous cities Uniform across centers (base)19,709Population Pro-rata population Uniform across centers (+3.5%)21,037Benchmarks Cases Pro-rata active cases Uniform across centers (+10.5%)23,622Locations Optimized: minimizes deaths Uniform across centers (+24.0%)25,615Full Optimized: minimizes deaths Optimized: minimizes deaths (+34.5%)Optimized: minimizes deaths Optimized: minimizes deaths 23,000Optimization Proposed Fair: inter-state equity Fair: inter-center equity (+20.8%)
A downside of these two solutions (“locations optimized” and “fully optimized”), however, isthat they can result in very inequitable outcomes between states and between vaccination sites.The last (“proposed”) solution circumvents this challenge by imposing all fairness constraints(Equations (25)–(27)), and choosing tight values for the corresponding hyperparameters Θ L , θ V and θ P . As such, the proposed solution ensures fairness across various dimensions (site locations, vaccineallocation across states, and vaccine allocation across sites). Remarkably, even when constrainingthe optimization as such, the resulting proposed solution still results in 20% benefits, as comparedto the top-cities benchmark—saving 4,000 deaths over the three-month horizon.In comparison, the benchmarks cannot reach the same impact of vaccinations. The top-citiesand population-based benchmarks achieve similar outcomes, with 19,000–20,000 lives saved ascompared to a no-vaccination baseline. The case-based baseline performs better, by saving 21,000lives (10% more than the top-cities baseline). Yet, these numbers remain significantly lower thanthose achieved with our solution. The top-cities and population-based baselines, by relying ondemographic information exclusively, fail account for disparities in the severity of the across thecountry. The case-based baseline captures the status of the pandemic, but not its dynamics, thustreating similarly a state where the pandemic is already waning and a state where it is risingfast—although vaccinations have a stronger impact in the latter state than the former.These results underscore the main takeaways from this paper: optimization provides significantbenefits, as compared to simple benchmarks based on readily-available information at the outset ofthe vaccination campaign. Instead, the optimized approaches design vaccine distribution strategies ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? based on both demographics and epidemiological dynamics. The edge of optimization can be sig-nificant, by saving an extra 20% of lives with the same vaccine capacity. Stated differently, underthe proposed optimization approach, each vaccine is effectively “worth” 1.2 vaccines.A natural question to ask is whether optimization induces strong geographic disparities, bymerely shifting the benefits of the vaccines from one state to another. To explore this question,Figure 6 plots the number of vaccines distributed to each state, as compared to its population share(Figure 6a) and the number of lives saved, as compared to the top-cities benchmark (Figure 6b). lll lll l ll l lllll ll l llll l ll llll ll l lll ll l llll ll ll l ll ll ll ll l ll l llll ll lll lll l l lll ll lll l ll ll l l llll l lll l l Less than share of populationMore than share of populationVaccines per state (a) Vaccines allocation lll lll l ll l lllll ll l llll l ll llll ll l lll ll l llll ll ll l ll ll ll ll l ll l llll ll lll lll l l lll ll lll l ll ll l l llll l lll l l
Less than baselineMore than baselineLives saved per state (vs. population−based baseline) (b) Lives saved
Figure 6 Vaccine allocation (vs. population share) and lives saved (vs. population-based baseline) per state.
Recall that the optimization imposes fairness constraints in vaccine allocation (Equations (25)–(27)); yet, the remaining flexibility can be used strategically to target the populations where vac-cines can have the strongest impact. In fact, as Figure 6a shows, vaccines do not get distributedproportionally to each state’s population. For instance, states like Texas, Florida and New Yorkreceive a higher share of vaccines, whereas states like California, Pennsylvania and Illinois receivea lower share. This is expected, given the significantly higher number of lives saved under theoptimized (“proposed”) solution as compared to the population-based benchmark (Table 2).However, Figure 6b shows that these disparities in vaccine allocation do not result in sharpdisparities in public health outcomes. Specifically, the optimized solution saves hundreds of extralives (as compared to the population-based benchmark) in seven states with very different epi-demiological profiles. Texas benefits the most from optimization (with an additional 1,450 livessaved), followed by New York (440), Florida (380) and Iowa (290). At the same time, the optimizedsolution does not increase the death toll in any state by more than 100. Pennsylvania is the most ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites? negatively impacted state, with an estimated 85 additional deaths—well within the margin of errorof DELPHI–V. Overall, the proposed optimization approach can thus distribute vaccine capacitiesto combat the pandemic in some critical states without hurting others. A core challenge in vaccine distribution lies in the significant uncertainty regarding the dynamicsof the pandemic and the effect of vaccinations. To address this challenge, we assess the sensitivityand the robustness of the optimized solution when the structure of the DELPHI–V model andsome of its key parameters are perturbed. For each perturbation, we compare three solutions:– the top-cities baseline, evaluated with the new perturbed model– the re-optimized solution, optimized and evaluated with the new perturbed model– the proposed solution, optimized with the original model and evaluated with the new per-turbed model. Specifically, we first run the optimization with the initial inputs. We then introduceperturbations, and re-optimize vaccine allocation decisions (variables C and V ), while fixing thevaccination sites (variable x ). Indeed, in practice, vaccination sites are determined once and forall, whereas vaccine allocation can be re-adjusted as information becomes available.We then compare the number of deaths under these three solutions, estimated with the perturbedDELPHI–V model. Comparisons between the top-cities benchmark and the re-optimized solutionestimate the sensitivity of the benefits of optimization to the perturbations. Comparisons betweenthe top-cities benchmark and the proposed solution estimate the robustness of the solution. % i n c r e a s e i n s a v e d li v e s proposed solution (effectiveness = 0.9)re-optimized solution (actual effectiveness) (a) Vaccine effectiveness % i n c r e a s e i n s a v e d li v e s proposed solution (budget = 1 million)re-optimized solution (actual budget) (b) Vaccine budget Figure 7 Sensitivity and robustness of results with varying vaccine effectiveness and vaccine budget.
We first vary the two main drivers of the vaccination campaign, that are the vaccine effectiveness(parameter β ) and the vaccine budget (parameters B t ). Figure 7 reports the percent-wise increasein saved lives, as compared to the top-cities benchmark, for both the proposed and the re-optimized ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? solutions. The main takeaways fall under three categories. First, the benefits of optimization in-crease with vaccine effectiveness. This is expected, as a higher vaccine effectiveness increases theimpact of all strategies (in the extreme example where β = 0, all strategies have the same nullperformance). Second, the benefits of optimization decrease with vaccine budget. Indeed, in theextreme scenario with an infinite vaccine budget, any distribution strategy can immediately endthe pandemic, leaving essentially no space for optimization. As the budget becomes more scarce,the decisions of who receives a vaccine and when become increasingly complex, so the optimizedstrategy has a positive and significant impact on the spread of the disease. It is interesting to notethat this monotonic relationship would get inverted in the other regime with a small vaccine bud-get (again, all solutions perform identically when the budget gets to zero). Yet, with the currentvaccine capacities, the benefits of optimization are very strong, saving an extra 15–35% of lives.The third takeaway from Figure 7 is the high degree of robustness of the proposed solution. Inall but one experiment, the proposed solution remains optimal under the perturbed parameters(indicated by exactly the same benefits obtained with the proposed solution and the re-optimizedsolution). In the last case (with a daily budget of 500,000 vaccines), the proposed solution isdominated by the re-optimized solution, but remains within 2% of the new optimum. Obviously,the new values of vaccine effectiveness and vaccine budgets impact downstream vaccine allocationsand the dynamics of the pandemic. However, the location of vaccination sites is highly robust tovariations in vaccine effectiveness and vaccine budgets.Next, we test the robustness of the proposed solution to the dynamics of the pandemic. Oneassumption of DELPHI–V is that the infection rate (captured by the parameters α j and γ jt ) isidentical across age groups. However, serological evidence suggests potential disparities in infectionsacross sub-populations. To test this, we vary the infection rates with the risk class k , by adjustingthe relative variations based on the serological estimates from the US Center for Disease Control(2020a). Another assumption is that vaccines prevent mortality but not infection and transmission.In practice, vaccines may still reduce the risk of infection and the propensity to transmit the disease.To test this, we perturb the DELPHI–V model by assuming that a fraction β of vaccinated peopletransition directly to the immune state M (as opposed to the susceptible and vaccinated state S (cid:48) ).Figure 8 reports the results from these new perturbations, using the same nomenclature asFigure 7. The main takeaway is identical, in that the proposed solution remains near-optimal underthe two perturbations. In both tests, the proposed solution results in less than 50 extra deathsthan the re-optimized solution—a very small amount in view of the 20,000 lives saved by theoptimization. In other words, the proposed solution is not only robust to vaccination characteristics(effectiveness and budget) but also to the dynamics of the COVID-19 pandemic. ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? uniform age-dependentinfection rate0510152025 % i n c r e a s e i n s a v e d li v e s proposed solutionre-optimized solution (a) Age-dependent infection rates transmission no transmissionCan vaccinated people transmit?05101520 % i n c r e a s e i n s a v e d li v e s proposed solutionre-optimized solution (b) No transmission Figure 8 Robustness of results with age-dependent infection rates and no transmission from vaccinated people.
Finally, we vary two key parameters in the DELPHI–V dynamics: the infection and mortalityrates. Although these parameters are fitted from historical data, there remains considerable uncer-tainty regarding the dynamics of the pandemic and people’s behaviors through a period of massvaccination. Therefore, we introduce a random perturbation in the infection or mortality rate, ineach state. Specifically, we define 50 perturbation scenarios; in each one, we sample each parameterin each state independently, following a uniform distribution centered around the nominal valueand spanning ± ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites? F r e q u e n c y (a) Infection rates F r e q u e n c y (b) Mortality rates Figure 9 Robustness of the benefits of optimization with infection and mortality rates. at the downstream level. This, again, illustrates the effects of the non-linear SEIR dynamics onthe spread of the disease, and how we can leverage an epidemiological model such as DELPHI–Vto allocate resources strategically in order to combat the pandemic most effectively.From a practical standpoint, the results from this section are highly significant. Indeed, theDELPHI–V model, like any epidemiological model, only provides a rough approximation of the dy-namics of the pandemic and the effects of vaccinations. Yet, the robustness of the solution providesguarantees that the locations of mass vaccination sites, even optimized against this approximatemodel, remain highly robust when evaluated against alternative dynamics. This is obviously notto say that we can commit to a full-scale vaccine distribution strategy that spans the entire vac-cination campaign. However, the “here-and-now” decisions (i.e., the location of vaccination sites)are likely to remain near-optimal (or even optimal) in the next phases of the pandemic, ultimatelyenabling the vaccination campaign to save an extra 15–35% of lives.
7. Conclusion
This paper has presented a new prescriptive approach to optimize vaccine distribution strategies inresponse to the COVID-19 pandemic. The approach starts with a state-of-the-art epidemiologicalmodel called DELPHI, which augments SEIR models by capturing dynamics specific to COVID-19(under-detection, governmental and societal response, and declining mortality rates). This paperhas first proposed an extension, named DELPHI–V, which captures the effects of vaccinations andreflects the disaggregated impact of COVID-19 on mortality across risk classes (e.g., age groups).Then, this paper has embedded the predictive DELPHI–V model into an optimization model,termed DELPHI–V–OPT, to support vaccine distribution. DELPHI–V–OPT is formulated as abilinear optimization model, and solved using a tailored algorithm based on coordinate descent. ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites? We applied the model and algorithm to one of the priorities of the new Biden administration inthe United States: determining the locations of mass vaccination sites across the country. We for-malize the problem by selecting locations strategically to minimize the death toll of the pandemic,subject to practical and fairness constraints. Experimental results using real-world data suggestthat the proposed optimization approach can yield significant benefits, as compared to benchmarksolutions that locate vaccination sites based on readily-available demographic and epidemiologicalinformation. Specifically, the model can increase the effectiveness of the vaccination campaign by20%, saving an extra 4,000 lives over a three-month period. Remarkably, the proposed solutionachieves critical fairness objectives—by significantly reducing the death toll of the pandemic inseveral states without hurting others—and is highly robust to uncertainties and forecast errors—byachieving similar benefits under a vast range of perturbations.Obviously, the optimization approach developed in this paper is not without limitations. Forinstance, our experiments have only partitioned the population according to age groups, thusignoring other objectives such as prioritizing allocations to healthcare workers, other essentialworkers, or patients with comorbidities. Moreover, our epidemiological model does not captureheterogeneity in population mixing across subpopulations (e.g., interactions may be more intensein urban areas and between young people than in rural areas and between older populations).Similarly, our model only captures the first-order effects of vaccinations, ignoring for instanceheterogeneity across multiple vaccine types and double-dose vaccines. Finally, our methodologyrelies on a time discretization approximation and a coordinate descent approach, which do notyield theoretical guarantees on solution quality.Whereas these limitations undoubtedly motivate further research, this paper lays the first data-driven brick on the optimal distribution of COVID-19 vaccines. At a time where vaccine devel-opment and vaccine production are going full speed, the results from this paper highlight thecritical role of vaccine distribution strategies to combat the pandemic. Obviously, it is essentialto make every effort possible to develop new vaccines, enhance vaccine effectiveness, and produceas many vaccines as possible. But this paper highlights another lever that can be pulled to curbthe effect of the pandemic: strategically managing vaccine stockpiles to prevent the spread of thepandemic and mitigate its impact. As such, this paper can provide critical decision-making supportto governmental agencies as they are currently planning mass vaccinations around the globe.
Acknowledgments
The authors gratefully acknowledge T. Greg McKelvey Jr. for feedback and discussions that havegreatly improved the quality of the paper. ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites? References
Aaby K, Herrmann JW, Jordan CS, Treadwell M, Wood K (2006) Montgomery county’s public healthservice uses operations research to plan emergency mass dispensing and vaccination clinics.
Interfaces
Health Care Management Science
Management Science
Pro-ceedings of the National Academy of Sciences arXiv preprintarXiv:2012.12263 .Bandi H, Bertsimas D (2020) Optimizing influenza vaccine composition: From predictions to prescriptions.
Machine Learning for Healthcare Conference , 121–142 (PMLR).Basta NE, Chao DL, Halloran ME, Matrajt L, Longini Jr IM (2009) Strategies for pandemic and sea-sonal influenza vaccination of schoolchildren in the United States.
American journal of epidemiology arXiv preprint arXiv:2006.16509 .Chick SE, Mamani H, Simchi-Levi D (2008) Supply chain coordination and influenza vaccination.
OperationsResearch
Manu-facturing & Service Operations Management
PLoS One
The Journal of Ambulatory Care Management
Manufacturing & Service Operations Management
European journal of epidemiology ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites?
European Journalof Operational Research
PLoS Med
American Journal of Epidemiology
Science
Management science
Nature nanotechnology
New England Journalof Medicine .Graham BS (2020) Rapid COVID-19 vaccine development.
Science
Scientific reports
New England journal of medicine .Jacobson SH, Sewell EC, Proano RA (2006a) An analysis of the pediatric vaccine supply shortage problem.
Health Care Management Science
Vaccine
Epidemics
Nature
Science .Kornish LJ, Keeney RL (2008) Repeated commit-or-defer decisions with a deadline: The influenza vaccinecomposition.
Operations Research ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites?
Bulletin of mathematical biology
Chemical Engineering Research and Design medRxiv
URL .Longini Jr IM, Ackerman E, Elveback LR (1978) An optimization model for influenza a epidemics.
Mathe-matical Biosciences
NewEngland Journal of Medicine
Management Science medRxiv .Matrajt L, Halloran ME, Longini Jr IM (2013) Optimal vaccine allocation for the early mitigation of pan-demic influenza.
PLoS Comput Biol
Science
JAMA
Journal of theoretical biology medRxiv .Rastegar M, Tavana M, Meraj A, Mina H (2021) An inventory-location optimization model for equitableinfluenza vaccine distribution in developing countries during the covid-19 pandemic.
Vaccine medRxiv . ertsimas, Digalakis Jr., Jacquillat, Li, Previero: Where to locate mass vaccination sites?
NatureNanotechnology
Operations Research
Mathematical biosciences
Service Science
OR spectrum
Journal of Applied Probability
Operations Research
European Journal of Operational Research ertsimas, Digalakis Jr., Jacquillat, Li, Previero:
Where to locate mass vaccination sites? Appendix