Optimizing Spatio-Temporal Allocation of the COVID-19 Vaccine Under Different Epidemiological Landscapes

An efficient and safe vaccine is expected to allow people to return to normal life as soon as possible. However, vaccines for new diseases are likely to be in short supply during the initial deployment due to narrow production capacity and logistics. There is an urgent need to optimize the allocation of limited vaccines to improve the population effectiveness of vaccination. Existing studies mostly address a single epidemiological landscape. The robustness of the effectiveness of other proposed strategies is difficult to guarantee under other landscapes. In this study, a novel vaccination allocation model based on spatio-temporal heterogeneity of epidemiological landscapes is proposed. This model was combined with optimization algorithms to determine the near-optimal spatio-temporal allocation for vaccines with different effectiveness and coverage. We fully simulated the epidemiological landscapes during vaccination, and then minimized objective functions independently under various epidemiological landscapes and degrees of viral transmission. We find that if all subregions are in the middle or late stages of the pandemic, the difference between the effectiveness of the near-optimal and pro-rata strategies is very small in most cases. In contrast, under other epidemiological landscapes, when minimizing deaths, the optimizer tends to allocate the remaining doses to sub-regions with relatively higher risk and expected coverage after covering the elderly. While to minimize symptomatic infections, allocating vaccines first to the higher-risk sub-regions is near-optimal. This means that the pro-rata allocation is a good option when the subregions are all in the middle to late stages of the pandemic. Moreover, we suggest that if all subregions are in the period of rapid virus transmission, vaccines should be administered to older adults in all subregions simultaneously, while when the epidemiological dynamics of the subregions are significantly different, priority can be given to older adults in subregions that are still in the early stages of the pandemic. After covering the elderly in the region, high-risk sub-regions can be prioritized.

Introduction

The rapid spread of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has triggered a public health and economic crisis worldwide. As of January 15, 2022, there have been more than 3.1 billion cases and 5.5 million deaths reported (1). To combat this crisis, authorities have implemented various non-pharmaceutical interventions according to local conditions. Both empirical analysis and mathematical modeling suggest that non-pharmacological interventions, particularly lockdown and strict quarantine measures (2–4), are essential to mitigate the spread of the virus in the short term. However, given the high socio-economic costs of such measures, a long-term solution—an effective and safe vaccine-remains urgently needed. With the successive implementation of mass vaccination campaigns, the huge demand for vaccines has left many countries and regions facing severe resource shortages and supply imbalances (5, 6). Optimizing the allocation of limited vaccines is an urgent and critical issue for all countries.

The mathematical models help to inform the public and policymakers about possible scenarios for the development of infectious diseases and the potential effectiveness of different intervention methods. The compartmental models can be easily used to simulate various interventions. Epidemiological parameters in these models and simulation results are easy to understand. Based on the characteristics of this COVID-19 pandemic, a series of extensions for the classical SIR model have been proposed, such as age-group separation (7, 8), the introduction of asymptomatic infections (9, 10), deaths (11), and immunity period (12), and the incorporation of graph-based spatial components (10, 13). The efficacy and effectiveness of most vaccines on various endpoints have been assessed through rigorous controlled clinical trials or observational studies, such as BNT162b2 (14, 15), ChAdOx1 (16), mRNA-1273 (17), and inactivated SARS-CoV-2 vaccine (18). The optimal allocation of vaccines involves three main dimensions: object, timing, and places. Most of such studies focused on to whom and suggested that vaccine should be prioritized to people by age and risk. For the individuals in a region, prioritized vaccination of the elderly is beneficial in reducing severe symptomatic infections and deaths, while vaccinating younger age groups first can minimize symptomatic infections (11, 12, 19–21). The research on the spatio-temporal allocation of vaccines can be divided into two parts. Some studies used mathematical models to evaluate the effectiveness of different allocation strategies (22–26). Another approach was using optimization algorithms to find the possible optimal strategy (27–29). In existing studies, only the allocation for the epidemiological landscape (mainly refers to epidemic dynamics in sub-regions in our work) of a region over a certain period was identified. This does help to address the vaccine allocation under the corresponding epidemiological landscape. However, epidemiological landscapes change over time, and different countries may face different epidemiological landscapes when deploying vaccines. Once the epidemiological landscape changes, the effectiveness of the proposed vaccine allocation strategy will be difficult to guarantee, i.e., there may be other better allocation strategies (see the Supplementary Materials). The effect of spatio-temporal heterogeneity of epidemiological landscapes on vaccine allocation is ignored to some extent in previous studies. These studies lack the optimal deployment of vaccines under different epidemiological landscapes.

The current vaccine spatial allocation strategy usually allocates available doses to sub-regions in proportion to the size of their population (referred to as the pro-rata allocation). For policymakers, this strategy may be easier to implement. However, such a strategy that ignores epidemiological landscapes may not always maximize the population effectiveness of vaccines. In this view, we proposed a new model of vaccination allocation using a graph-based spatial model and an extended SIR-based temporal model. We simulated the epidemic dynamics in subregions during vaccine deployment sufficiently through this model to set up a variety of possible vaccination scenarios. The optimization algorithm was then run under different scenarios to independently minimize two metrics of infection and disease burden: proportion of cumulative infections and deaths prevented. We also evaluated vaccine allocation in combination with non-pharmacological interventions. In epidemiology, a non-pharmaceutical intervention (NPI) is any method to reduce the spread of an epidemic disease without requiring pharmaceutical drug treatments. Vaccination is a pharmaceutical intervention by definition. However, we believe that vaccination is often not administered alone, but in conjunction with other non-pharmacological interventions. The non-pharmaceutical interventions may affect the effectiveness of vaccination strategies. There, we introduced vaccination at different intensities of non-pharmaceutical interventions to investigate the variation patterns of optimal strategies in the context of different degrees of viral transmission.

Methods

Vaccination Allocation Model

The vaccine allocation model consists of two components: an age-structured deterministic compartmental model to describe the temporal dynamics of infectious diseases within the sub-region and a graph-based spatial model to simulate the epidemiological landscapes that the region may face at the time of vaccination.

Age-Structured Deterministic Compartmental Model

We stratified individuals by 10-year age groups, in line with previous parameter estimation and data. We assumed subjects started susceptible to infection (S) and could become exposed but not yet infectious (E) after effective contact. After a latent period, exposed individuals developed an asymptomatic infection (A) or a symptomatic infection (I). Both symptomatic and asymptomatic infectious individuals had a certain probability of being detected (ID). Infectious individuals eventually recovered (R) or died (D) depending on the severity of their symptoms. Among them, asymptomatic individuals did not die but recovered at a given rate. We also assumed that recovered individuals became susceptible again at a given rate, reflecting eventual loss of temporary immunity from the infection (30, 31). We used the previously estimated (32) contact matrix C with age structure and corrected for reciprocity (33) based on the existing demographic structure of the subregion (34), adapting it to each subregion.

We introduced a leaked vaccine into the model and tracked the vaccinated individuals by subscript v. Vaccination was available to susceptible and recovered individuals. We assumed that the vaccine can have two effects on the vaccinated individuals. First, they are less likely to be infected with the virus compared to unvaccinated individuals (this effect is expressed by V_e). Second, the vaccine can also reduce the probability of developing symptoms upon infection (referred to as V_p). Vaccine effect on COVID-19 disease is defined as reduction in the likelihood to develop symptomatic diseases upon exposure (referred to as V_dis). The relationship (35) among V_dis, V_e, and V_p is shown as follows:

$\begin{array}{ll}{V}_{dis}=1-((1-V_e)(1-V_p))& (1)\end{array}$

A flowchart of the model is presented in Figure 1. The equations used in the model are as follows:

$\begin{array}{ll}{c}_{ij}=w_h{c}_{ij|home}+w_w{c}_{ij|work}+w_s{c}_{ij|school}+w_o{c}_{ij|other}& (2)\end{array}$

$\begin{array}{ll}{\lambda }_i=\beta m_i{\displaystyle \sum _{j=1}^9}{c}_{ij}(\partial (A_j+{\text{A}}_{v,j})+I_j+I_{v,j})/N_j& (3)\end{array}$

$\begin{array}{ll}\frac{d{S}_i}{dt}=-{\lambda }_i{S}_i+\mu (R_i+R_{v,i})& (4)\end{array}$

$\begin{array}{ll}\frac{d{E}_i}{dt}={\lambda }_i{S}_i-\sigma E_i& (5)\end{array}$

$\begin{array}{ll}\frac{d{A}_i}{dt}=(1-p_i)\sigma E_i-({\gamma }_A+{\omega }_A)A_i& (6)\end{array}$

$\begin{array}{ll}\frac{d{I}_i}{dt}=p_i\sigma E_i-({\gamma }_I+{\omega }_I)I_i& (7)\end{array}$

$\begin{array}{ll}\frac{dI{D}_i}{dt}={\omega }_A{A}_i+{\omega }_I{I}_i-{\gamma }_{ID}I{D}_i& (8)\end{array}$

$\begin{array}{ll}\frac{d{R}_i}{dt}={\gamma }_A{A}_i+(1-{\delta }_i)({\gamma }_I{I}_i+{\gamma }_{ID}I{D}_i)-\mu R_i& (9)\end{array}$

$\begin{array}{ll}\frac{d{D}_i}{dt}={\delta }_i({\gamma }_I{I}_i+{\gamma }_{ID}I{D}_i)& (10)\end{array}$

$\begin{array}{ll}\frac{d{S}_{v,i}}{dt}=-(1-V_e){\lambda }_i{S}_{v,i}& (11)\end{array}$

$\begin{array}{ll}\frac{d{E}_{v,i}}{dt}=(1-V_e){\lambda }_i{S}_{v,i}-\sigma E_{v,i}& (12)\end{array}$

$\begin{array}{ll}\frac{d{A}_{v,i}}{dt}=(1-(1-V_p)p_i)\sigma E_{v,i}-({\gamma }_A+{\omega }_A)A_{v,i}& (13)\end{array}$

$\begin{array}{ll}\frac{d{I}_{v,i}}{dt}=(1-V_p)p_i\sigma E_{v,i}-({\gamma }_I+{\omega }_I)I_{v,i}& (14)\end{array}$

$\begin{array}{ll}\frac{dI{D}_{v,i}}{dt}={\omega }_A{A}_{v,i}+{\omega }_I{I}_{v,i}-{\gamma }_{ID}I{D}_{v,i}& (15)\end{array}$

$\begin{array}{ll}\frac{d{R}_{v,i}}{dt}={\gamma }_A{A}_{v,i}+(1-{\delta }_i)({\gamma }_I{I}_{v,i}+{\gamma }_{ID}I{D}_{v,i})-\mu R_{v,i}& (16)\end{array}$

$\begin{array}{ll}\frac{d{D}_{v,i}}{dt}={\delta }_i({\gamma }_I{I}_{v,i}+{\gamma }_{ID}I{D}_{v,i})& (17)\end{array}$

where c_ij is the number of individuals in age group j contacted by an individual in age group i, c_ij|home, c_ij|work, c_ij|school, and c_ij|other are the number of contacts at home, work, school, and other locations, respectively, and w_h, w_w, w_s, and w_o are the weight of the contact matrix of the four locations, respectively. λ_i is the force of infection for an individual in age group i, β is the transmission coefficient, μ_i is the relative susceptibility to infection for age-i, and α is the relative infectiousness of asymptomatic infections. 1/μ is the average length of immunity. 1/σ is the mean duration of latent period. p_i is the proportion of the symptomatic in age group i. 1/γ_A, 1/γ_I, and 1/γ_ID are the average recover time of the asymptomatic, symptomatic and detected, respectively. ω_A and ω_I are the detection rate (via contact tracing and testing) for asymptomatic and symptomatic infections, respectively. δ_i is the fatality of the disease for age-i individuals.

FIGURE 1

Figure 1. The flowchart of the transmission model.

Graph-Based Spatial Model

Due to the spatio-temporal heterogeneity and uncertainty in the spread of infectious diseases and the implementation of interventions, the epidemiological dynamics of sub-regions when deploying vaccines can be either similar (same early, middle or late stage of a pandemic) or very different (sub-regions at different times of the pandemic). We introduced into the model the times at which vaccination begins in sub-regions to set up a series of possible vaccination scenarios based on this phenomenon. Different scenarios represent different epidemiological landscapes, as shown in Figure 2. For sub-region k under scenario s, if vaccination starts on day $t_k^s$, it means that the virus has been transmitted for ($t_k^s$ – 1) days. Depending on the scenario, the times were randomly generated within a reasonable range. This will affect the initial status of each compartment, such as sub-regions that are in the early stage of the pandemic with a greater proportion of susceptible individuals. We quantified the intensity of the intervention using the control reproduction number (R_c, defined as the average number of secondary COVID-19 infections produced by typically infected individuals in a susceptible population with control measures). The R_c values for sub-regions were also randomly generated within a certain range. We used these R_c values to estimate some of the parameters in the model that were relevant to non-pharmacological interventions. This means that these parameters may be different for different sub-regions. Full simulation details can be found in subsection Vaccination Simulation Scenarios.

FIGURE 2

Figure 2. Example of the simulated prevalence of sub-regional symptomatic infections under different scenarios without vaccination. (A) Scenario 1 (All subregions are in the early stages of the pandemic). (B) Scenario 2 (All subregions are in the middle stages of the pandemic). (C) Scenario 3 (All subregions are in the late stages of the pandemic). (D) Scenario 4 (The epidemic dynamics of the subregions differ significantly).

We considered a baseline basic reproductive number R₀ = 2.5 (no interventions) to calculate the transmission coefficient β. We defined the next-generation matrix (NGM) as

$\begin{array}{ll}NG{M}_{ij}=\beta m_i{\displaystyle \sum _j}{c}_{ij}(\partial d_A+d_I)& (18)\end{array}$

R₀ is the spectral radius of the NGM (36). d_Aand d_I denote respectively the length of time individuals spend in states A and I.

Risk Assessment of Viral Transmission in Subregions

The effective reproduction number (R_e, defined as the average number of secondary cases per infectious case in a population made up of both susceptible and non-susceptible hosts) was usually used to measure epidemic transmission. Its value will change as the pandemic proceeds, as some people may have gained immunity through infection or vaccination. It is not very desirable to use the effective reproduction number at a given moment in isolation when prioritizing resources. Policymakers need to look at the whole picture and use it alongside other indicators. As shown in Equation 19, we assessed the risk of transmission in subregions by the effective reproduction number and the prevalence of symptomatic infections.

$\begin{array}{ll}{\text{r}}_k=R{S}_k+\eta I{S}_k& (19)\end{array}$

where r_k is the risk of viral transmission in the sub-region k. RS_k and IS_k are the total effective reproduction number (sum of effective reproduction number over T (our time horizon) days) and the total prevalence of symptomatic infections (sum of symptomatic infections over T days) without vaccination in the sub-region k, respectively, both of which are normalized. Considering that R_e could better reflect the transmission risk, we assumed η = 0.5.

Optimization Algorithm

The main focus of our study was on time and space. Within the subregion, our optimizer would prioritize vaccine allocation to the elderly. This is in line with most current policies and is more ethical. Two objective functions were ultimately chosen: the proportion of cumulative symptomatic infections and deaths prevented 1 year (our time horizon) after vaccination initiation compared with the unvaccinated base case. The reason why we choose symptomatic infections is that these individuals are more infectious compared to asymptomatic infections, and minimizing symptomatic infections is better for controlling the further spread of the virus. Moreover, symptomatic infections are closely related to public health policy concern: the burden on the health system. Therefore, we believe that although policies tailored to reduce deaths are of great significance, there is still a need to optimize the allocation of vaccines for symptomatic infections. This can provide more insight and flexibility in the formulation of public health policy.

$\begin{array}{ll}{\displaystyle \underset{X}{min}}f(X)& (20)\end{array}$

$\begin{array}{ll}s.t.\{\begin{array}{c}{\displaystyle \sum _k}{X}_k=1\\ X_{k}N\le P_k\\ 0\le X_k\le 1\end{array}& (21)\end{array}$

where f (X) is the objective function and X is the decision variable. X_k is the proportion of vaccines allocated to subregion k to the total available vaccines. P_k is the number of people in sub-region k and N is the total number of people in the region.

To increase the chances of approaching the global optimal solution, we combined coarse global search with genetic algorithm to explore the entire space of possible combinations of vaccine allocation. Before running the genetic algorithm, we performed coarse-grained grid search over the entire decision variable space, with a search step of 0.1. For example, a point in the grid is X = (0,0.1,0,0,0.1,0.4,0.2,0.2,0,0), which represents the sub-regions 2, 5, 6, 7, 8 are allocated to 10, 10, 40, 20, and 20% of the available vaccines, respectively. Each grid point needed to satisfy the above constraints, if not, the grid point would be discarded. We evaluated the objective function on all feasible grid points and selected the best 25 points. The best 25 points obtained above and the pro-rata allocation vector were put into the initial population of the genetic algorithm, and the remaining individuals of the initial population were randomly generated by the genetic algorithm according to the constraints. It is worth noting that under high vaccine supply (60−100%), the selection of feasible grid points may be <25 due to the reduction in the number of decision variables that meet the constraints. If the vaccine supply reaches 100%, the only theoretically feasible decision variable is the pro-rata vector. In this case, all individuals of the initial population were randomly generated by the algorithm itself.

Results

We partitioned the continental US into the 10 Standard Federal Regions (10 HHS Regions) established by the US Office of Management and Budget. The 10 HHS Regions are our study subjects. More information on the 10 HHS regions can be found in Supplementary Materials. Given the current V_dis data and the uncertainty in the specific V_e and V_p values, we evaluated the optimal use of 50 vaccines with specific effectiveness and coverage (V_dis ranging from 60 to 100% and vaccination coverage ranging from 10 to 100% of the total population, in in-crements of 10%). We set V_e = 50% while varying the value of V_p according to V_dis.

Our study included four categories of scenarios and three intensities of non-pharmacological interventions. We assumed that the intensities of the non-pharmacological interventions in subregions were close and set, [2.0,2.3], [1.5,1.8], and [1.2,1.4] for low, moderate, and high intensity non-pharmacological interventions, respectively. We assumed that sub-regions deployed vaccines at a certain rate per day (rollout speeds of 0.05% of their population per day) until the supply was exhausted. Since each country will have different vaccination rates, the at-once allocation would make the results more general. Therefore, we also used this method to determine the near-optimal strategy under different scenarios (see the Supplementary Materials). All simulations were performed on the Matlab 2021a platform. For the genetic algorithm, we used the global optimization toolbox provided by it. The implementation details and hyperparameters of the genetic algorithm are shown in Table 1.

TABLE 1

Table 1. The values for the hyperparameters of the GA.

Simulation Setup

Parameter Estimation

The values of most parameters in the model referred to previous studies, while some of the parameters related to non-pharmacological interventions were estimated using genetic algorithm. The decision variables of the genetic algorithm are all parameters to be estimated. The fitness value is the absolute value of the difference between the target R_c and the estimated R_c. The values of relevant parameters are shown in Table 2.

TABLE 2

Table 2. The values and source of parameters in the model.

Vaccination Simulation Scenarios

The possible scenarios were divided into two main categories: epidemic dynamics in sub-regions are relatively close or different significantly. For the first category of scenarios, we built three sub-scenarios (Scenarios 1, 2, and 3). Scenarios 1, 2 and 3 represent epidemiological landscapes with all subregions in the early, middle, and late stages of the pandemic, respectively, while scenario 4 represents that the epidemiological dynamics of the sub-regions are significantly different. The generation intervals of the times when sub-regions start vaccination under different scenarios are as follows:

$\begin{array}{llll}{t}_k^s& \in & \{\begin{array}{l}(p_a{t}_k^p,(p_a+l_w)t_k^p),s=1\\ ((1-\frac{1}{2}{l}_w)t_k^p,(1+\frac{1}{2}{l}_w)t_k^p),s=2\\ ((1+p_a)t_k^p,(1+p_a+l_w)t_k^p),s=3\\ (0,t_k^l),s=4\end{array}& (22)\end{array}$

$\begin{array}{ll}{l}_w=\{\begin{array}{ll}0.1,& s=1,2,3\\ 1,& s=4\end{array}& (23)\end{array}$

where $t_k^s$ represents the time when vaccination begins in sub-region k under scenario s. $t_k^p$ represents the expected timing of the pandemic peak in sub-region k. l_w is the scaling factor of the generation interval to control the degree of proximity of the epidemic dynamics among sub-regions. p_a (0 < p_a < 1–l_w) is used to determine the location of the interval. In scenarios 1 and 3, we adjusted the value of p_a and conducted multiple experiments to check whether the near-optimal vaccine strategy changed when the interval changed. $t_k^l$ is the pandemic duration for sub-region k.

Uncertainty and Sensitivity Analysis of Model Parameters

We assumed that each parameter obeys a uniform distribution over its range of values. Latin hypercube sampling (LHS) was performed for the entire parameter space. The model was run 1,000 times to analyze the sensitivity of response function values to changes in a single input parameter using the partial rank correlation coefficient (PRCC) as an indicator. The results are shown in Figure 3.

FIGURE 3

Figure 3. Sensitivity analysis of the main parameters of the model. Each plot shows the sensitivity of the corresponding response function value to changes in a single main input parameter of the model. The green (red) bars indicate sensitive (insensitive) and gray bars indicate results that are not statistically significant.

These results suggest that COVID-19 can be controlled effectively in the population by some interventions such as wearing masks, maintaining social distance, or large-scale vaccination. We chose model parameters that were more sensitive to the objective functions as parameters to be estimated: weights of contacts at home, school and other locations, and detection rates for asymptomatic as well as symptomatic infections.

Near-Optimal Vaccine Allocation Differs for Different Epidemiological Landscapes

The near-optimal spatio-temporal allocation for the same scenario is nearly identical under different intensities of non-pharmacological interventions. In this section, we analyzed how the near-optimal strategy to minimize symptomatic infections changed for different epidemiological landscapes and present results under low-intensity non-pharmacological interventions.

Compared with no vaccination, both strategies substantially mitigated viral transmission in all scenarios, especially when vaccine coverage is relatively high. If V_dis= 60% and vaccine coverage was 50%, the near-optimal strategy could avert 40% (UI: 37−43%), 17% (UI: 16−17%), 42% (UI: 33−57%), and 52% (UI: 43−62%) of symptomatic infections in the four scenarios on average, respectively (Figure 4). This suggests that mass vaccination can go a long way in alleviating the spread of COVID-19 and also emphasizes the importance of distributing vaccines as soon as possible. Because the population effectiveness of vaccination will diminish as the pandemic proceeds. The near-optimal strategy under scenarios 1 and 4 could provide more gain compared to other scenarios, averting up to 9% more (V_dis = 60%, ~60% of vaccine coverage) and 15% more symptomatic infections (V_dis = 60%, ~40% of vaccine coverage) at the greatest difference between the two strategies, respectively. In scenario 1, when the vaccine coverage was ~30−70%, it was optimal to vaccinate higher-risk sub-regions at high coverage, while the pro-rata allocation was close to the near-optimal strategy for other coverages (Figure 5A). In stark contrast, in Scenarios 2 and 3, the pro-rata strategy was close to optimal regardless of vaccine coverage and effectiveness (Figure 5B). The near-optimal vaccine allocation is similar for scenarios 4 and 1, with the difference that the priority was also given to higher-risk sub-regions when vaccine coverage is less than 30% (Figure 5C). This may be since that in scenario 1, the risk of transmission in the sub-regions is not very different and the virus is rapidly spreading. All subregions require large amounts of vaccines to slow the further spread of the virus. While in scenario 4, some sub-regions are already in the middle and late stages of the pandemic. In the case of insufficient vaccines, the available doses can be allocated to the sub-regions in the early stages of the pandemic. Both strategies tended to perform similarly as vaccine effectiveness and coverage increased.

FIGURE 4

Figure 4. Percentage of cumulative symptomatic infections averted for both strategies under different scenarios. When V_dis = 60% and low-intensity non-pharmacological intervention, the percentage of cumulative symptomatic infections averted for the near-optimal strategy (green) and the pro-rata strategy (purple). Bars represent the mean of multiple experimental results and error bars represent uncertainty intervals (UI).

FIGURE 5

Figure 5. Near-optimal allocation strategies to minimize symptomatic infections under scenario 1 (A), 3 (B) and 4 (C). For each heat map, each row from left to right is the decreasing direction of transmission risk, representing the total vaccine supply (percentage of the total population vaccinated) and each column represents a different subregion. Colors represent the percentage of the population in a sub-region to be vaccinated.

Near-Optimal Vaccine Allocation Changes With Objective Functions

Next, we investigated the effect of the vaccination objectives on the near-optimal strategy and present some of the results under moderate intensity non-pharmacological interventions. The near-optimal vaccine allocation for the two objective functions differed more in scenarios 1 and 4. While in scenarios 2 and 3, whatever the objective function was, the pro-rata allocation was close to optimal in most cases.

In Scenario 1, for low vaccine coverage, the near-optimal strategy to minimize symptomatic infections prioritized the coverage in higher-risk sub-regions; while when minimizing deaths, the optimizer's allocation results were similar to the pro-rata strategy. As more vaccines became available (30−60% vaccine coverage), the near-optimal strategy for symptomatic infections remained the same, while the near-optimal strategy for minimizing deaths first covered the elderly in all sub-regions as much as possible and then allocated the remaining dose to the subregion with the expected higher risk and coverage (Figures 6A,C). The main reason for this phenomenon is that the main contributors to the spread of the virus in this COVID-19 pandemic are the younger age groups, but the mortality rate in the elderly is much higher than in other age groups. In our optimizer, the vaccine is first distributed to the elderly in the sub-region. Allocating vaccine to older people in the region first helps to minimize deaths. The effective control of the epidemic requires mass vaccination among young people. The near-optimal spatio-temporal strategy under scenario 4 is overall similar to that of scenario 1, except that with fewer vaccines available, the optimization algorithm may only focus on allocating vaccines to sub-regions in the early and middle stages of the pandemic (Figures 6B,D). In addition, we found that the gain from optimizing vaccine allocation to improve pro-rata strategy was relatively small when minimizing the deaths (Figure 7). This may be because prioritizing older adults is more beneficial for minimizing deaths within a sub-region than other vaccine prioritization strategies.

FIGURE 6

Figure 6. Near-optimal allocation strategies for different objective functions. Here, we only show the results for V_dis = 60 and 80%. (A) The near-optimal strategies to minimize symptomatic infections in scenario 1. (B) The near-optimal strategies to minimize symptomatic infections in scenario 4. (C) The near-optimal strategies to minimize deaths in scenario 1. (D) The near-optimal strategies to minimize deaths in scenario 4. For each heat map, each row from left to right is the decreasing direction of transmission risk, representing the total vaccine supply (percentage of the total population vaccinated) and each column represents a different subregion. Colors represent the percentage of the population in a sub-region to be vaccinated.

FIGURE 7

Figure 7. Percentage of cumulative symptomatic infections (A) and deaths (B) averted for both strategies. Solid lines with solid circles and dashed lines with hollow circles represent the near-optimal and pro-rata strategies, respectively. For clarity, we highlight 1 set of results in Figure 6 in each scenario. The shaded areas represent the range of effectiveness of strategies under multiple epidemiological landscapes simulated.

Effect of the Intensity of Non-pharmacological Interventions on Vaccine Allocation

Since the epidemic dynamics of the subregions in each experiment were randomly generated, the epidemiological landscapes under scenario 4 may vary considerably under different intervention intensities. Therefore, we only discuss the effect of the intensity of non-pharmacological interventions on the near-optimal vaccine allocation for the first 3 scenarios.

We found that the intensity of non-pharmacological interventions had a greater effect on the near-optimal strategy under scenario 1. As the intensity of non-pharmacological interventions increased, the effectiveness of both strategies improved to varying degrees compared with no vaccination (Figure 8). For example, when V_dis = 60% and the total vaccine coverage was 50%, the near-optimal strategy in scenario 1 could prevent 43% (UI: 37−43%) of symptomatic infections if low-intensity non-pharmacological interventions were implemented along with vaccination. The symptomatic infections averted under moderate or high intensity non-pharmacological interventions would increase to 76% (UI: 51−76%) and 99% (UI: 96−99%), respectively (Figure 8A). While the effectiveness of the pro-rata strategy increased from 36% (UI: 34−36%) to 71% (UI: 45−71%) and 98% (UI: 87−99%; Figure 8A). The implementation of high-intensity non-pharmaceutical interventions not only increases the population effectiveness of vaccination but also requires fewer vaccines to control the epidemic. Vaccinating everyone in scenario 1 would still be difficult to completely control the outbreak under low-intensity non-pharmaceutical interventions. However, under high-intensity interventions, only 50% of the population would need to be near-optimally vaccinated to control the pandemic. Furthermore, when non-pharmacological interventions were intensified, the optimizer tended to allocate vaccines more evenly across all sub-regions. For the same vaccine effectiveness and coverage, the range of optimal vaccine coverage for sub-regions decreased with increasing intensity of the non-pharmacological interventions, especially when the goal of vaccination was to minimize cumulative symptomatic infections and vaccine effectiveness and coverage were not sufficiently high (Figure 9).

FIGURE 8

Figure 8. Percentage of cumulative symptomatic infections (A) and deaths (B) prevented under different intensities of interventions. When V_dis = 60% and vaccine coverage is 10−100% of the population, the percentage of cumulative symptomatic infections and deaths prevented for the near-optimal (solid line with solid circles) and pro-rata (dashed line with hollow circles) strategy. For clarity, we highlight 1 set of results in each scenario. The shaded areas represent the range of effectiveness of strategies under multiple epidemiological landscapes simulated.

FIGURE 9

Figure 9. Near-optimal vaccine coverage in sub-regions. Each violin represents the distribution of the primary optimal vaccine coverage for all sub-regions in scenario 1 when minimizing the number of symptomatic infections for a certain intensity of non-pharmaceutical intervention.

Discussion

The COVID-19 pandemic is not only a serious threat to human life and health but also has a huge impact on economic development and social stability. Non-pharmacological interventions are effective in reducing the incidence of COVID-19 cases. However, in the long term, these measures are difficult to implement consistently and large-scale vaccination is highly desirable. However, during the actual deployment of the vaccine, the demand for the vaccine may far outweigh the supply. Therefore, optimizing the allocation of limited vaccines makes good practical sense and can provide valuable information for policy formulation. In this work, we proposed a novel vaccine spatio-temporal allocation model that can reproduce the epidemiological landscapes that the region may face when deploying vaccines. We used this model and optimization algorithms to determine near-optimal vaccine spatiotemporal allocation strategy in different scenarios. Our findings suggest that the performance of vaccine spatiotemporal allocation strategies is driven by several complex factors, including vaccine effectiveness and coverage, epidemic dynamics of subregions when allocating, and non-pharmacological interventions implemented.

Based on the experimental results, we provide the following recommendations for vaccination policy. If all subregions are in the middle or late stage of the pandemic when the vaccine is deployed, the existing pro-rata strategy is desirable, regardless of vaccine effectiveness and coverage. This is because under such epidemiological landscapes, for both vaccination targets, the gains from optimizing vaccine allocation are minimal compared to the pro-rata allocation. When sub-regions are in the same period of rapid virus transmission or when the epidemic dynamics of sub-regions differ significantly, prioritizing the higher-risk sub-regions may be more beneficial in controlling further virus transmission. While to minimize the deaths in the region, the elderly in all subregions need to be vaccinated as soon as possible. Under scenarios 1 and 4, the gain from optimizing vaccine allocation to improve pro-rata strategy was relatively small when minimizing the deaths. When the vaccine coverage is less than 30%, the effectiveness of the two strategies is almost the same under scenario 1. Based on this phenomenon, we suggest that older adults in all subregions should be vaccinated simultaneously if all subregions are in a period of rapid virus transmission, while older adults in subregions still in the early stages of a pandemic could be vaccinated first when the epidemiologic dynamics of each subregion differ significantly. After covering older adults in that region, priority can be given to high-risk subregions. In addition, our results suggest that the synergistic implementation of multiple pandemic interventions can lead to better pandemic control and that enhanced non-pharmaceutical interventions can reduce vaccination pressure. If non-pharmaceutical interventions are enhanced when deploying vaccines, the optimization strategy will tend to allocate the vaccine more evenly across subregions. This means that vaccines can be allocated more evenly across sub-regions based on the above recommendations to further enhance the effectiveness of vaccination if public health officials can encourage or urge the public to maintain good self-protection habits while vaccinating.

It is worth noting that this work still has certain limitations. First, the model assumed that the vaccine provided equal protection to all vaccinated individuals and that the effect was constant over the simulated time. Nevertheless, there may be some differences in the protective effect of the vaccine in different age groups, and asymptomatic infected individuals may result in weaker protection (39). The protective effect of artificial immunity induced by vaccination may decline over time (40). If immunity is transient, then these results will apply only for that duration. Second, the same rollout speed may be less desirable, so different rollout speeds for subregions may be the focus of our subsequent study, such as a faster rollout of vaccines in sub-regions with more supply.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://zenodo.org/record/6574113#.You7Ke5BxPY.

Author Contributions

WC and JZ designed research and conceived the experiments. WC, JZ, and HD conducted the experiments and analyzed the results. WC, JZ, and XW contributed to the drafting of the work. WC, JZ, XW, XT, YT, and ZM contributed to the review and editing of the manuscript. All authors read and approved the final manuscript.

Funding

This study was funded by the National Key Research and Development Program of China (2018YFB0505304) and the National Natural Science Foundation of China (Grant No. 41671409).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpubh.2022.921855/full#supplementary-material

Abbreviations

COVID-19, coronavirus disease 2019; NPIs, non-pharmaceutical interventions; R₀, basic reproduction number; R_e, effective reproduction number; R_c, control reproduction number; NGM, next-generation matrix.

References

1. Johns Hopkins University and Medicine. Coronavirus COVID-19 Global Cases by the Center for Systems Science and Engineering (CSSE) at Johns Hopkins University (JHU). (2019). Available online at: https://coronavirus.jhu.edu/map.html (accessed July 16, 2021).

2. Anderson SC, Edwards AM, Yerlanov M, Mulberry N, Stockdale JE, Iyaniwura SA, et al. Quantifying the impact of COVID-19 control measures using a Bayesian model of physical distancing. PLoS Comput Biol. (2020) 16:e1008274. doi: 10.1371/journal.pcbi.1008274

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Cao W, Dai H, Zhu J, Tian Y, Peng F. Analysis and evaluation of non-pharmaceutical interventions on prevention and control of COVID-19: a case study of Wuhan City. ISPRS Int J Geo-Inf. (2021) 10:480. doi: 10.3390/ijgi10070480

CrossRef Full Text | Google Scholar

4. Peak CM, Kahn R, Grad YH, Childs LM Li R, Lipsitch M, et al. Individual quarantine versus active monitoring of contacts for the mitigation of COVID-19: a modelling study. Lancet Infect Dis. (2020) 20:1025–33. doi: 10.1016/S1473-3099(20)30361-3

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Khamsi R. If a coronavirus vaccine arrives, can the world make enough. Nature. (2020) 580:578–80. doi: 10.1038/d41586-020-01063-8

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Usher A D. COVID-19 vaccines for all? Lancet. (2020) 395:1822–3. doi: 10.1016/S0140-6736(20)31354-4

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Lazebnik T, Bunimovich-Mendrazitsky S. The signature features of COVID-19 pandemic in a hybrid mathematical model—implications for optimal work–school lockdown policy. Adv Theory Simul. (2021) 4:2000298. doi: 10.1002/adts.202000298

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Bunimovich-Mendrazitsky S, Stone L. Modeling polio as a disease of development. J Theor Biol. (2005) 237:302–15. doi: 10.1016/j.jtbi.2005.04.017

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Lazebnik T, Bunimovich-Mendrazitsky S, Shami L. Pandemic management by a spatio–temporal mathematical model. Int J Nonlinear Sci Numer Simul. (2021). doi: 10.1515/ijnsns-2021-0063

CrossRef Full Text | Google Scholar

10. Viguerie A, Lorenzo G, Auricchio F, Baroli D, Hughes TJ, Patton A, et al. Simulating the spread of COVID-19 via a spatially- resolved susceptible–exposed–infected–recovered–deceased (SEIRD) model with heterogeneous diffusion. Appl Math Lett. (2021) 111:106617. doi: 10.1016/j.aml.2020.106617

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Bubar KM, Reinholt K, Kissler SM, Lipsitch M, Cobey S, Grad YH, et al. Model-informed COVID-19 vaccine prioritizationstrategies by age and serostatus. Science. (2021) 371:916–21. doi: 10.1126/science.abe6959

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Foy BH, Wahl B, Mehta K, Shet A, Menon GI, Britto C. Comparing COVID-19 vaccine allocation strategies in India: a mathematical modelling study. Int J Infect Dis. (2021) 103:431–8. doi: 10.1016/j.ijid.2020.12.075

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Lazebnik T, Alexi A. Comparison of pandemic intervention policies in several building types using heterogeneous population model. Commun Nonlinear Sci Numer Simul. (2022) 107:106176. doi: 10.1016/j.cnsns.2021.106176

CrossRef Full Text | Google Scholar

14. Dagan N, Barda N, Kepten E, Miron O, Perchik S, Katz MA, et al. BNT162b2 mRNA Covid-19 vaccine in a nationwide mass vaccination setting. N Engl J Med. (2021) 384:1412–23. doi: 10.1056/NEJMoa2101765

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Tartof SY, Slezak JM, Fischer H, Hong V, Ackerson BK, Ranasinghe ON, et al. Effectiveness of mRNA BNT162b2 COVID-19 vaccine up to 6 months in a large integrated health system in the USA: a retrospective cohort study. Lancet. (2021) 398:1407–16. doi: 10.1016/S0140-6736(21)02183-8

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Voysey M, Clemens SA, Madhi SA, Weckx LY, Folegatti PM, Aley PK, et al. Single-dose administration and the influence of the timing of the booster dose on immunogenicity and efficacy of ChAdOx1 nCoV-19 (AZD1222) vaccine: a pooled analysis of four randomised trials. Lancet. (2021) 397:881–91.

PubMed Abstract | Google Scholar

17. Baden LR, El Sahly HM, Essink B, Kotloff K, Frey S, Novak R, et al. Efficacy and safety of the mRNA-1273 SARS-CoV-2 vaccine. N Engl J Med. (2021) 384:403–16. doi: 10.1056/NEJMoa2035389

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Jara A, Undurraga E A, González C, et al. Effectiveness of an inactivated SARS-CoV-2 vaccine in Chile. N Engl J Med. (2021) 385:875–84. doi: 10.1056/NEJMoa2107715

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Jentsch PC, Anand M, Bauch CT. Prioritising COVID-19 vaccination in changing social and epidemiological landscapes: a mathematical modelling study. Lancet Infect Dis. (2021) 21:1097–106. doi: 10.1016/S1473-3099(21)00057-8

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Matrajt L, Eaton J, Leung T, Dimitrov D, Schiffer JT, Swan DA, et al. Optimizing vaccine allocation for COVID-19 vaccines shows the potential role of single-dose vaccination. Nat Commun. (2021) 12:1–18. doi: 10.1038/s41467-021-23761-1

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Matrajt L, Eaton J, Leung T, Brown ER. Vaccine optimization for COVID-19: who to vaccinate first? Sci Adv. (2021) 7:abf1374. doi: 10.1126/sciadv.abf1374

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Araz OM, Galvani A, Meyers LA. Geographic prioritization of distributing pandemic influenza vaccines. Health Care Manag Sci. (2012) 15:175–87. doi: 10.1007/s10729-012-9199-6

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Azman A S, Lessler J. Reactive vaccination in the presence of disease hotspots. Proc R Soc B Biol Sci. (2015) 282:20141341. doi: 10.1098/rspb.2014.1341

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Wu JT, Riley S, Leung GM. Spatial considerations for the allocation of pre-pandemic influenza vaccination in the United States. Proc R Soc B Biol Sci. (2007) 274:2811–7. doi: 10.1098/rspb.2007.0893

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Krzysztofowicz S, Osińska-Skotak K. The use of GIS technology to optimize COVID-19 vaccine distribution: a case study of the city of Warsaw, Poland. Int J Environ Res Public Health. (2021) 18:5636. doi: 10.3390/ijerph18115636

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Zhou S, Zhou S, Zheng Z, Lu J. Optimizing spatial allocation of COVID-19 vaccine by agent-based spatiotemporal simulations. GeoHealth. (2021) 5:e2021GH000427. doi: 10.1029/2021GH000427

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Venkatramanan S, Chen J, Fadikar A, Gupta S, Higdon D, Lewis B, et al. Optimizing spatial allocation of seasonal influenza vaccine under temporal constraints. PLoS Comput Biol. (2019) 15:e1007111. doi: 10.1371/journal.pcbi.1007111

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Molla J, Chávez AP, Hiraoka T, Ala-Nissila T, Kivelä M, Leskelä L. Adaptive and optimized COVID-19 vaccination strategies across geographical regions and age groups. PLoS Comput Biol. (2022) 18:e1009974. doi: 10.1371/journal.pcbi.1009974

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Lemaitre JC, Pasetto D, Zanon M, Bertuzzo E, Mari L, Miccoli S, et al. Optimal control of the spatial allocation of COVID-19 vaccines: Italy as a case study. medRxiv. (2022). doi: 10.1101/2021.05.06.21256732

CrossRef Full Text | Google Scholar

30. Kojima N, Klausner J D. Protective immunity after recovery from SARS-CoV-2 infection. Lancet Infect Dis. (2022) 22:12–4. doi: 10.1016/S1473-3099(21)00676-9

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Sariol A, Perlman S. Lessons for COVID-19 immunity from other coronavirus infections. Immunity. (2020) 53:248–63. doi: 10.1016/j.immuni.2020.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Prem K, Cook AR, Jit M. Projecting social contact matrices in 152 countries using contactsurveys and demographic data. PLoS Comput Biol. (2017) 13:e1005697. doi: 10.1371/journal.pcbi.1005697

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Arregui S, Aleta A, Sanz J, Moreno Y. Projecting social contact matrices to different demographic structures. PLoS Comput Biol. (2018) 14:e1006638. doi: 10.1371/journal.pcbi.1006638

PubMed Abstract | CrossRef Full Text | Google Scholar

34. CDC, Population. Projections. Available online at: https://wonder.cdc.gov/wonder/help/populationprojections.html (accessed June 26, 2021).

Google Scholar

35. Mehrotra DV, Janes HE, Fleming TR, Annunziato PW, Neuzil KM, Carpp LN, et al. Clinical endpoints for evaluating efficacy in COVID-19 vaccine trials. Ann Int Med. (2021) 174:221–8. doi: 10.7326/M20-6169

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. doi: 10.1016/S0025-5564(02)00108-6

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Davies NG, Klepac P, Liu Y, Prem K, Jit M, Eggo RM. Age-dependent effects in the transmission and control of COVID-19 epidemics. Nat Med. (2020) 26:1205–11. doi: 10.1038/s41591-020-0962-9

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Laxminarayan R, Wahl B, Dudala SR, Gopal K, Mohan BC, Neelima S, et al. Epidemiology and transmission dynamics of COVID-19 in two Indian states. Science. (2020) 370:691–7. doi: 10.1126/science.abd7672

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Long QX, Tang XJ, Shi QL Li Q, Deng HJ, Yuan J, et al. Clinical and immunological assessment of asymptomatic SARS-CoV-2 infections. Nat Med. (2020) 26:1200–4. doi: 10.1038/s41591-020-0965-6

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Khoury DS, Cromer D, Reynaldi A, Schlub TE, Wheatley AK, Juno JA, et al. Neutralizing antibody levels are highly predictive of immune protection from symptomatic SARS-CoV-2 infection. Nat Med. (2021) 27:1205–11. doi: 10.1038/s41591-021-01377-8

PubMed Abstract | CrossRef Full Text | Google Scholar