Introduction

Supply- or deployment-limited SARS-CoV-2 vaccines [1] pose the urgent question of a fair distribution of the available doses [2]. Current prioritization approaches typically target groups at higher risk of severe outcomes [3, 4], or their indirect protection by vaccinating those with higher disease transmission [3, 5, 6]. In this paper we will discuss how taking into account spatial heterogeneity in disease transmission when designing prioritization strategies significantly improves the effectiveness of vaccination campaigns. The distribution of doses inside each country is limited by the logistic capabilities of the healthcare network and the rate at which the vaccine stock is replenished. Decisions concerning the best allocation strategies are to be taken under these constraints. Moreover, both the complex coupling between regions due to human mobility and the spatial heterogeneity in disease history and control interventions make the discovery of such optimal allocation strategies an arduous task.

We propose an optimal control framework to explore COVID-19 vaccine distribution in space and time. We study the SARS-CoV-2 epidemic in Italy, where strong spatial effects arise from the geography of the disease, heterogeneous lockdown exit strategies, and post-lockdown control measures [7]. The optimal control framework is applied to a spatial model that has proved its reliability for Italy [8, 9], whose parameters are here sequentially updated through the assimilation of a year-long epidemiological record. This allows us to unravel the best possible vaccination strategy and probe the impact of vaccine allocations over the 107 Italian provinces.

The problem of vaccine allocation is of primary importance for public-health officials, epidemiologists, and economists [10, 11]. Roll-out strategies are conventionally based on the prioritization of individuals at risk, such as health workers and elderly people [12–15]. However, the heterogeneous ways in which different regions may be affected by each successive wave raise questions about spatial prioritization strategies. What is the best feasible spatial allocation, given supply and logistic constraints? Would that differ significantly from current non-spatially-optimized plans? Should vaccines be distributed based on demography or would it be better to prioritize areas currently subject to an outbreak? How relevant are the susceptibility profile and modeled future transmission in each region?

Epidemiological modeling has long been used to answer questions about the impact of vaccination campaigns, often by comparing outcomes under different scenarios [16, 17]. Optimization of epidemiological models, i.e, the search for the best possible course of action that maximizes or minimizes an objective metric, has been carried out theoretically since the seventies [18–20]. Recent dramatic improvements of both algorithms [21] and computational power prompted applied studies using different methods to rigorously find optimal mitigation strategies [22–24]: most of the time trough iterative parameter search [25, 26], but also using genetic algorithms [27], greedy algorithms [28] or solving the Hamilton-Jacobi-Bellman equations [29, 30].

Interesting developments have recently arose during the ongoing SARS-CoV-2 pandemic [13, 31, 32]. The urgency of effective vaccination campaigns led to the development of modeling frameworks for the optimization of vaccine allocation, based on age or risk [3, 4, 12, 13], space [33], dose timing [34, 35], and the deployment of testing resources, using optimal control [36] or Bayesian experimental design [37], along with prioritization based on social contact networks [38].

To the best of our knowledge, optimal spatial allocation of COVID-19 vaccines at a country scale has never been performed yet. This question is distinct from, and complementary to, risk-based prioritization. Spatial heterogeneities in disease transmission are complex, as seen during the initial outbreaks [8, 9, 39], supporting the significance of the posed problem towards an effective control of the epidemic. However, the connectivity network underlying spatial epidemiological models may generate complex large-scale control problems whose solution requires tailored formulations and efficient algorithms.

This work aims to find optimal strategies for this problem through state-of-the-art optimization methods based on distributed direct multiple shooting, automatic differentiation, and large-scale nonlinear programming [40–43]. This approach allows us to solve the large-scale optimization problems arising from epidemiological models, even when considering hundreds of spatial nodes.

Materials and methods

The formulation of the optimal control problem has three main components: 1) an objective function to be minimized, here the number of new infections; 2) the spatial epidemiological model [8, 9] governing the transmission dynamics with the daily vaccination rates in each province as control variables; and 3) the set of constraints that the control must satisfy, in our case the limitations on vaccine administration rate in each province and the total vaccine stock in Italy.

1) Objective function. Optimizing calls for a metric, whose selection is critical in determining the optimal solution and its outcome. The choice of an objective function relates to health, economy, and ethics. Possible candidates are the minimization of, e.g., DALYs (the Disability-Adjusted Life Years), the number of deaths, disease exposure, or economic loss [44]. All these objectives are linked and may be combined. As the model considered for this work does not have risk classes, we optimize for the minimization of the incident infections in Italy from January 4, 2021, to April 4, 2021. Minimization of the deaths would yield the same results under the assumptions used in the model.

2) Epidemiological model. Incidence and deaths are projected using the spatially distributed epidemiological model devised by Gatto et al. [8] and further improved by Bertuzzo et al. [9]. The model subdivides the Italian population into 107 provinces represented as a network of connected nodes. Each province has local dynamics describing the number of individuals present in each of the model compartments: susceptible S, exposed E, pre-symptomatic P (incubating infectious), symptomatic infectious I, asymptomatic infectious A, hospitalized H, quarantined Q, recovered R, and dead D. A tenth compartment, vaccinated individuals V, is added to the original nine, as shown in Fig 1A. Compartments P, A, and I have different degrees of infectiousness and contribute to the force of infection (Eq (C) and Eq (D) in S1 Text), which represents the rate at which susceptibles S become infected and enter the exposed compartment E.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Model setup.

(A) Diagram representing the compartments of the epidemiological model and the possible transitions in a single province. We control the vaccination rate (teal arrows), aiming at minimizing incident infections (pink arrow). Individuals in compartments outside of the yellow block are able to move along the mobility network shown in (B), hence the force of infection in a province is coupled with the dynamics of other connected provinces. To reduce the problem to a tractable size, we only consider the most important connections (red edges) when optimizing, but we use the full network (red and grey edges) to assess our strategies. A discussion on the effect of this simplification is provided in section C.2 in S1 Text. Nodes’ size and color display each province’s population, and edges’ width shows the strength of the coupling between a pair of provinces. Base map layer from the Istituto Nazionale di Statistica; Istat, istat.it, CC-BY 3.0.

https://doi.org/10.1371/journal.pcbi.1010237.g001

Except for those in H, Q, D, or I states, a fraction of individuals commutes between provinces along the mobility network, thus we introduce node-to-node disease transmission along the network shown in Fig 1B. Thus, the force of infection in each province has a local and a mobile component. The local component describes transmission among the individuals that do not leave the province. The mobile component considers that local susceptibles may enter in contact with infected individuals that are traveling, and oppositely, susceptible commuters may become infected through contact with local infected. Connected provinces contribute to this process depending on the strength of the mobility fluxes from and to the node of interest. These mobility fluxes change in time due to the governmental policies introduced to reduce transmission among regions (complete model equations and more details about the data used to construct the mobility network and its use in the model are presented section B in S1 Text).

The epidemiological model, previously calibrated during the first wave of COVID-19 in Italy [8, 9], is updated up to January 4, 2021 using an iterative particle filtering, which infers the regional transmission on a moving temporal window of two weeks. This data assimilation scheme allows us to capture the second wave of infections that hit Italy in the Fall of 2020, a necessary requirement to generate model projections that take into account the whole epidemic history, as shown in Fig 2. In our approach, model projections are described by an ensemble of a hundred trajectories associated with different parameters, whose distributions quantify the model uncertainty. We consider two projection scenarios characterized by two possible rates of epidemic transmission, see Fig 2. The “Optimistic” scenario assumes a constant lowering of transmission from January 4, 2021, to April 4, 2021; the “Pessimistic” scenario considers a gradual increase in transmission until mid-February 2021, which results in a third wave.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Data assimilation and scenarios for optimization.

Comparison between the model outputs (95% confidence interval (CI) of the ensemble, blue shaded area) and the corresponding epidemiological data (red circles, obtained from the bulletins of the Dipartimento della Protezione Civile, https://github.com/pcm-dpc/COVID-19) from March 2020 to January 2021. The orange and green shaded areas respectively show the ensemble dynamics (95% CI) of what we called pessimistic and optimistic transmission scenarios from January to April. The optimal vaccination strategy in the optimistic (or pessimistic) scenario is computed for the the continuous green (or orange) line, representing the model trajectory obtained using the median of each ensemble parameter. (A) The data on the daily hospitalizations is estimated as described in [9]; this data at the regional level is assimilated on a moving window of 14 days to update the model parameters describing the local transmission rates (see section D in S1 Text). (B) Daily number of newly exposed individuals versus the reported positive cases. Note that the large discrepancy between model and data during the first wave is due to the low testing capacity at the beginning of the epidemic (C) Daily number of deaths.

https://doi.org/10.1371/journal.pcbi.1010237.g002

The control variable is the vaccination rate in each province. We assume one-dose vaccines with an instantaneous 100% efficacy, while in reality the vaccine efficacy and immunity duration depend on the vaccine type. As we focus on spatial patterns and differences among vaccination strategies for a given supply, this assumption does not affect the conclusion of our work. Finally, we impose that vaccine protection persists during the three months of projection considered.

For each scenario, the optimal control problem is solved for one reference model trajectory, whose parameters and state on January 4, 2021, are obtained as the median values of the 100 model realizations. In this way, the reference trajectory approximately represents the ensemble median in each province. Then, we assess the effectiveness of the optimal allocation on the full ensemble of trajectories.

3) Constraints. We define two types of constraints: supply constraints, which determine the weekly delivery to the national stockpile; and logistic constraints, which limit the maximum rate of local vaccine distribution in each province.

The supply constraints ensure that the model does not distribute more vaccines than what is actually available in stock. We assume that the national supply of vaccine doses is empty on January 4, 2021, and is replenished every Monday. We consider four scenarios with weekly deliveries of 125’00, 250’000, 479’700 (the latter is the most realistic and our baseline value), and 1M vaccine doses (additional results for scenarios with 1.5M and 2M doses delivered each week are shown in Table B and Fig I in S1 Text).

From the national stockpile, doses may be allocated to any of the 107 Italian provinces, but the logistic constraints limit the rate at which it is possible to distribute the vaccines in each province. We assume that the maximum number of individuals who can be vaccinated in a province per day is proportional to the province’s population, such that the national maximum distribution capacity equals 500’000 doses per day, i.e., 3.5M per week if every province vaccinates at its maximum rate (which in retrospect is close to Italy’s vaccination rate as of May 1, 2021).

The objective, the model, and the constraints may be tailored to specific applications within the proposed framework.

Using state-of-the-art large-scale nonlinear optimization solvers and automatic differentiation, we solve each scenario (optimistic and pessimistic, with different weekly stockpile deliveries) for the optimal vaccines allocation.

Results

We obtain the optimal vaccination strategies for a set of eight scenarios drawn from the spatial model from January 4, 2021, to April 4, 2021. These scenarios are a combination of two projection scenarios (pessimistic vs optimistic) and four assumptions on the weekly stockpile delivery (125’000, 250’000, 479’700, or 1M doses delivered per week). In each scenario, the optimal solution is a spatially explicit vaccine roll-out policy, i.e. an indication of the number of vaccine doses to be deployed in each province each day.

Performance of the optimal and alternatives vaccination strategies

Spatial prioritization based on epidemiological criteria, such as past [16] or future [17] incidence, has often been used in both real campaigns and prospective studies. In order to measure the improvements yielded by the optimal allocation strategy, we compare it against 12 alternative approaches which distribute the available weekly vaccine doses among provinces. These alternative strategies use an indicator variable to rank provinces, either i) their population; ii) the number of susceptible individuals (per inhabitant or absolute) at the beginning of the projection; iii) the future incidence as projected by the epidemiological model (per inhabitant or absolute) or iv) constant, equal for all provinces. The incidence indicator rankings are updated every day to reflect the change caused by past decisions. For each indicator we propose two variants: after ranking all the provinces, strategies either focus on the provinces where the indicator is the largest or allocate to all provinces proportionally to the indicator. Additionally, we further consider the greedy strategy presented in [28]. In the main text we present the results of the optimal strategy, the second-best strategy overall (indicator: incidence per inhabitant; allocation: focused), and proportional allocation strategies for incidence, susceptibility, and population. All other results, along with detailed pseudo-code for each strategy are presented in section E and F in S1 Text.

For each of the eight scenarios considered, we compute the number of averted infections with respect to a zero-vaccination baseline, and the number of averted infections per vaccination dose (see Table 1). In the optimistic transmission scenario, characterized by a recess of the epidemic, the vaccination campaign has a lower impact on the averted infections per dose as only a small percentage of the vaccinated individuals would have been at risk of transmission. As expected, the impact of the vaccination campaign is more evident in the pessimistic scenario where the optimal strategy averts up to 2.54 million infections given weekly stockpile deliveries of one million doses. By virtue of the law of diminishing returns, the number of averted infections per dose decreases (from 0.413 to 0.196) when increasing the weekly stockpile.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Absolute number of averted infections (in millions) and averted infections per dose during the first three months of 2021 as evaluated for the reference trajectory (see Fig 2) for each strategy.

The first column represents the considered scenarios of weekly stockpile replenishment, i.e., the number of doses delivered to Italy every week, ranging from 125’000 to one million.

https://doi.org/10.1371/journal.pcbi.1010237.t001

The optimal solution always outperforms all the explored alternative strategies in terms of the number of averted infections and in terms of averted infections per dose allocated (see Table 1, and Tables A and B in S1 Text). After extensive hyper-parameter tuning, the alternative strategy focusing on the provinces with the largest incidence has the closest results to the optimal strategy, with a difference of less than 10% in each scenario. Instead, the other strategies are significantly less effective. The improvement between optimal and incidence-based (proportional) allocation is significant, ranging from 9.0% (pessimistic, 1M doses/week) to 27.4% (optimistic, 125’000 doses/week). In Fig 3, the black diamonds represent the percentage of averted infections obtained using each strategy for the reference trajectory, normalized with respect to the averted infections resulting from the optimal strategy. We observe that, in both the optimistic and pessimistic scenarios, the optimal strategy has the largest relative benefits for the smallest stockpile.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Comparison between different vaccine allocation strategies.

Percentage of averted infections per vaccine dose from January 4, 2021, to April 4, 2021, resulting from province-scale vaccine allocation strategies for both the pessimistic (A) and the optimistic (B) scenarios based on the following vaccination strategies: the optimal solution, proportional to the province population, proportional to the susceptible individuals, proportional to the projected incidence, and focused on the provinces with the largest weekly incidence (see color codes in the legend). We optimize the vaccine allocation for the reference trajectory (the median trajectory in the model projections, indicated as diamonds in the figure), and assess the performance of the computed vaccination strategy over the whole posterior of trajectories (boxen plots). For each projection scenario, the results are normalized by the number of averted infections in the reference solution (see Table 1 for the absolute values). Results for alternative scenarios and vaccination strategies are shown in Fig I in S1 Text.

https://doi.org/10.1371/journal.pcbi.1010237.g003

In the pessimistic scenario (see Fig 3A), when 479’700 doses are available each week, the averted infections associated with the optimal strategy in the reference projection are 0.272 per dose: 24.6% more compared to the strategies based on population or susceptible distributions (0.205 averted infections per dose), and more than 14% higher compared to the strategy based on the projected incidence per inhabitant (0.232 averted infections per dose), while only 4% higher than the focused incidence strategy. These differences are smaller but still significant when increasing the weekly stockpile deliveries up to 1M doses; similar results are obtained also for the optimistic transmission scenario (Fig 3B).

We recall that the optimal control strategy is computed for a reference model trajectory, which is the median of an ensemble of 100 realizations. To further investigate the effectiveness of the optimal solution, we apply it to all trajectories of the ensemble. The box plots in Fig 3 display the main quantiles of the averted infections computed for the ensemble of trajectories. We observe that the optimal strategy still yields better results on the ensemble of projections compared to the other strategies, thus suggesting that the computed solution is robust even under the presence of perturbations in the forecasts of the epidemic dynamics. More importantly, for each realization of the ensemble and each projection scenario, the optimal strategy systematically averts more infections than any of the other control strategies.

The same conclusions hold while considering all alternative strategies (see Fig I and Table A and B in S1 Text). Finally, we present in Fig O in S1 Text, an additional sensitivity analysis where the epidemiological dynamics are shuffled between provinces, and thus are completely different from the projections used to design the strategies. This confirms that, despite its specialization, the optimal strategy does not under-perform with respect to alternatives in such scenarios.

Our results therefore suggest that it is possible to considerably increase the impact of vaccination campaigns by optimizing the vaccine allocation in space and time. For this task, optimal control provides the best possible strategy and sets a benchmark for the theoretical potential of a vaccination campaign.

Analysis of the optimal vaccine allocation

The optimal vaccine allocation obtained as the solution of the optimal control framework is complex to analyze. We will attempt to unravel the mechanism behind its performance. The strategy must obey the imposed logistic and supply constraints: 1) The vaccine stockpile is replenished every Monday by a fixed amount of doses (e.g., 479’700 doses in the baseline scenario), and 2) the maximum possible distribution capacity per province is limited, proportionally to the province population, so that the number of doses distributed across the country can be of 0.5M per day at maximum (more details in Fig A in S1 Text).

We display the optimal vaccination strategy in time for 479’700 doses/week in the pessimistic scenario in Fig 4. We observe that the optimal allocation respects the two constraints on distribution (Fig 4A) and supply (Fig 4B). We observe that no province is vaccinated at the maximum possible rate during the whole campaign and that provinces display a variety of vaccination patterns. We also note that all vaccines received every Monday are always fully distributed during the following week, but that the rate of delivery on a national level increases with time (Fig 4B). Surprisingly, the optimal solution favors precise targeting over the speed of delivery, in order to allocate more doses to provinces where the impact of vaccines on the whole system is projected to be higher. Hence, in order to control infections, precise targeting may trump delivery speed.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Optimal vaccine allocation for the baseline, pessimistic scenario.

(A) Cumulative proportion of vaccine doses administered in the 107 provinces, some of which are highlighted. The local distribution rate is limited by a rate that is proportional to the population. This logistic constraint is visualized here as the maximum slope, equal for every province. (B) Stacked cumulative absolute number of vaccines in the 107 provinces of Italy. The national stockpile is shown in black and is replenished every week (on Mondays) with 479’700 doses. We display the name of the provinces with a final allocation of more than 150’000 doses.

https://doi.org/10.1371/journal.pcbi.1010237.g004

Furthermore, we observe in the optimal solution that every time a province is vaccinated, the rate of vaccination is equal to the maximum rate allowed by the local logistic constraint, as it is the case for any focused alternative vaccination strategies. In Fig 5, one can already see by visual inspection that the optimal allocation distributes most of the available doses to a few provinces with high incidence. These provinces are neither the most connected nor the most populous in Italy. The optimal strategy makes then use of the information on the network connectivity to fine-tune the allocation and spreads the doses across more provinces than the incidence-based strategy.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Spatial distribution patterns for the optimal allocation (left) and alternative strategies based on population, incidence, and susceptibility (additional alternative strategies are presented in section E and F in S1 Text).

We show, for each province and strategy, the proportion of vaccinated individuals after the implementation of the strategy (top), the number of averted infections per inhabitant with respect to a no vaccination baseline (middle), and the proportion of individuals who are still susceptible at the end of the control horizon (bottom) for the pessimistic transmission scenario with a weekly stockpile delivery of 479’700 doses. Base map layer from the Istituto Nazionale di Statistica; Istat, istat.it, CC-BY 3.0.

https://doi.org/10.1371/journal.pcbi.1010237.g005

To further investigate these patterns, in Fig 6A we display the number of administered doses versus the incidence projected without vaccines (the proxy variables leading to the second-best control performance), both normalized according to the resident population in each province. We observe an allocation pattern whereby provinces with a higher incidence receive more vaccines, for both the pessimistic and optimistic scenarios. However, the allocation is nonlinear with respect to the projected incidence, suggesting that in order to better control the epidemic the optimal allocation strategy takes into account other factors such as the importance of each province within the mobility network, as well as the proportion of susceptibles. This complexity is further highlighted in the heatmaps in Fig 6B, where a complete picture of the spatio-temporal patterns is unraveled. For a stockpile delivery of 125’000 doses per week, we observe different patterns of local allocation in the provinces that receive vaccines. Along the course of the epidemic, the optimal strategy varies which groups of provinces are most vaccinated. When the weekly stockpile delivery is increased, we observe that the optimal solution allocates new doses by both further reinforcing already prioritized provinces (acting as a focused strategy) and by vaccinating new provinces (acting as a proportional strategy). The ability to allocate each dose where it is the most efficient, without tuning any hyper-parameter (such as a threshold when to act proportional versus focused) is a major benefit of the optimal control strategy. Additional analyses are provided in Figs J–N in S1 Text.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Analysis of the optimal solution.

(A) Vaccinated population according to the optimal strategy against the projected incidence without vaccination, both normalized by province population and considering the scenario with a weekly stockpile delivery of 479’700 doses. (B) Heatmap of optimal allocation in space (y-axis, one square per province in alphabetical order) and time (x-axis, one square per day) for increasing weekly delivery scenarios (left to right). The color represents the proportion of individuals vaccinated on this day in this province by the optimal solution for the pessimistic transmission scenario; black is the maximum logistic rate per inhabitant (which is equal for all provinces). The same graphs for the optimistic transmission scenario and other weekly deliveries are shown in Figs L and M in S1 Text.

https://doi.org/10.1371/journal.pcbi.1010237.g006

Discussion and conclusions

Without any supply constraint, each country would vaccinate its population as fast as possible according to the available infrastructure. However, limitations in vaccine supply and rate of delivery are a reality for every country, hence the available doses should be deployed in space and time following a fair and effective strategy.

In stockpile-limited settings, like most current vaccination campaigns worldwide, careful allocation may significantly increase the number of averted infections and deaths. The goal is to distribute the vaccines where they have the strongest beneficial impact on the dynamics of the epidemic. However, designing an algorithm capable of computing spatially optimal allocation strategies in real heterogeneous settings is far from trivial and our approach is, to the best of our knowledge, the first attempt in this direction.

We developed a novel optimal control framework that delivers the best vaccination strategy under realistic supply and logistics constraints. This allows us to compute the allocation strategy that maximizes the number of averted infections during a projection of the COVID-19 epidemic in Italy from January 4, 2021, to April 4, 2021. Our results show that the optimal strategy has a complex structure that mainly reflects the projected incidence of each province, but also takes into account the spatial connectivity provided by the mobility network and the landscape of acquired population immunity. Although the reason why this strategy is optimal is not immediately intuitive, our simulations clearly outline that it significantly outperforms other, more straightforward strategies. This comparison suggests that the simplicity underlying intuitive vaccination strategies may undermine their effectiveness, and calls for complementing these simple approaches with rigorous and objective mathematical tools, such as optimal control, that fully account for the complexity of the problem.

With the present work, we showed that it is possible to solve optimal control problems for spatially explicit dynamical models of infectious diseases at a national scale, thus overcoming the computational limitations that, up to now, precluded this kind of applications. The proposed framework can account for any compartmental epidemic model, with up to hundreds of connected spatial nodes. Supply and logistic constraints can be adapted to the actual landscape of decisions faced by the stakeholders, such as no/reduced vaccine delivery on weekends, or the need for fairness in vaccine distribution, e.g., by ensuring that each province receives at least a fixed fraction of the available vaccines. This is especially important as in the optimal allocation some provinces might receive no vaccine at all. Moreover, while we assumed single-dose vaccines, one could optimize the timing between doses of multi-dose vaccines, in addition to every other control, in the same framework.

The present work is not devoid of limitations. The main one is that the optimal vaccination strategy strongly depends on the projection of the underlying epidemiological model. These projections, as every epidemiological projection, are subject to several sources of uncertainty, especially for long horizons, e.g., due to model design and calibration [47], assumptions about future events in transmission scenarios, and unforeseen events that may change the course of the epidemic (such as the importation of cases, the emergence of new virus variants, changes in disease awareness or social distancing policies). These aspects have an impact on the optimal vaccination strategy, which is reliable only if the projections given by the underlying model dynamics are sufficiently accurate. A successful approach developed by the automatic control community to tackle that issue, named Model Predictive Control [48], consists in compensating for the performance losses expected over long horizons by constantly adapting the optimal strategy. In this context, Model Predictive Control can be implemented using the following steps: (a) at the beginning of each week, the state of the system is estimated by using newly acquired epidemiological data; (b) the optimization problem is solved over a fixed prediction horizon using the estimated state as an initial condition; (c) the optimal strategy for the first week is applied and, as soon as the next week starts, these steps are repeated starting from (a). This method corrects the model inaccuracies by constantly resetting the initial state to the estimated one. Additionally, constraints may be updated to account for unexpected deliveries or new orders. Future work will aim at further evaluating the benefits of implementing this scheme for the design of optimal vaccination strategies. Stochastic control is another possible direction of improvement: the presented method solves the optimal control problem for the median trajectory of an ensemble of projections, but it is theoretically possible to compute an optimal vaccine allocation that accounts for the whole uncertainty range—either optimizing for the best expected performance or the best worst-case performance. However, methods for robust or stochastic control, that entail feeding the whole posterior distribution of trajectories into the optimal control framework, are most likely not computationally tractable for our problem. Instead, we provide the evaluation across the whole posterior of the optimal solution for the median trajectory, and find that its performance remains superior to the one of alternative strategies. We observe that for our problem, this method is enough to outperform other strategies. Moreover, the sensitivity analysis provided in section H in S1 Text demonstrates that despite its specialization, the optimal allocation strategy does not perform comparatively poorly when the model projections are inaccurate.

The epidemiological model underlying our control optimization has known validity and limitations [8, 9]. A significant simplification is the use of Google Community Mobility Reports (see section D in S1 Text) to estimate the variations in mobility across provinces and as a proxy for changes in social contacts. This approach does not explicitly account for complex social interactions and interventions (such as opening/closure of schools, shopping, mask-wearing, and social distancing), and data assimilation is necessary to adjust modeled transmission to the observed one. An additional limitation of the model for the specific scopes of this work is that it does not account explicitly for risk-based classes, and thus does not account for the heterogeneity that may result from the demography of the population, as well as from the age-related transmission and clinical characteristics associated with COVID-19. While surely limiting for operational use of the tools, we note that the scope of this paper is to provide a proof of concept of the relevance of spatial effects, which have not been addressed so far in the literature. To that end, we are confident that our results support the relevance of the research question posed. Our framework can be extended to optimize across both spatial and risk heterogeneity, provided that sufficient computational capacities are available to solve the resulting optimal control problem.

A counter-factual assumption in this work is that we consider a one-dose vaccine with full and instantaneous efficacy against transmission. At the time of development, the details about COVID-19 vaccines were not released, and this hypothesis allowed us to demonstrate our framework in a simple setting. Our framework can be further extended to account also for the simultaneous deployment of different vaccine types, some of which may require the administration of two doses. This extension is the subject of ongoing research, in particular to extend the modeling tools described here to accommodate the peculiarities of each authorized vaccine candidate while designing effective spatio-temporal deployment strategies.

In conclusion, in this work we presented and analyzed the optimal allocations of vaccines against SARS-CoV-2 at a country scale under different scenarios of epidemic transmission and vaccine availability. We designed a general optimal control pipeline that performs large-scale nonlinear optimization of epidemiological controls. We demonstrated our framework by coupling it with an existing model of COVID-19 transmission over the 107 provinces of Italy. This model was updated with a data assimilation scheme to reflect the epidemic state and history as of January 2021. Within our optimal control framework, we discretized, transformed and simplified the model to find the best possible allocation of vaccines under realistic stockpiles and logistic constraints. The optimal solution outperforms other strategies by a significant margin and proves robust across the uncertainty of the underlying model. As such, besides inherent limitations, it provides a benchmark against which other, possibly simpler vaccine rollout strategies can be usefully compared. We analyzed the mechanisms behind our optimal allocation and concluded that the complex interplay of spatial heterogeneity and human mobility requires non-trivial prioritization strategies in order to achieve the maximum effectiveness.

Supporting information

Acknowledgments

JCL acknowledges Nicolas Richart for help on scaling up computations on the EPFL high-performance computing clusters.