How can age-based vaccine allocation strategies be optimized? A multi-objective optimization framework

Human life is deeply influenced by infectious diseases. A vaccine, when available, is one of the most effective ways of controlling the spread of an epidemic. However, vaccine shortage and uncertain vaccine effectiveness in the early stage of vaccine production make vaccine allocation a critical issue. To tackle this issue, we propose a multi-objective framework to optimize the vaccine allocation strategy among different age groups during an epidemic under vaccine shortage in this study. Minimizing total disease onsets and total severe cases are the two objectives of this vaccine allocation optimization problem, and the multistage feature of vaccine allocation are considered in the framework. An improved Strength Pareto Evolutionary Algorithm (SPEA2) is used to solve the optimization problem. To evaluate the two objectives under different strategies, a deterministic age-stratified extended SEIR model is developed. In the proposed framework, different combinations of vaccine effectiveness and vaccine production capacity are investigated, and it is identified that for COVID-19 the optimal strategy is highly related to vaccine-related parameters. When the vaccine effectiveness is low, allocating most of vaccines to 0–19 age group or 65+ age group is a better choice under a low production capacity, while allocating most of vaccines to 20–49 age group or 50–64 age group is a better choice under a relatively high production capacity. When the vaccine effectiveness is high, a better strategy is to allocate vaccines to 65+ age group under a low production capacity, while to allocate vaccines to 20–49 age group under a relatively high production capacity.

Introduction

Human development has been shaped by infectious diseases, and some infectious disease outbreaks have even led to the demise of a country, such as the biblical pharaonic plagues that hit Ancient Egypt in the middle of the Bronze Age at ~1715 BC, the plague in Athens from 430 to 425 BC that ended the golden age of Athens, and the Black Death bubonic plague in Europe in 1348, which killed millions of people (1). In recent decades, we have seen that SARS in 2003, H1N1 influenza in 2009, MERS in 2013, Ebola in 2014 and the COVID-19 since 2019 have led to millions of deaths, leading to substantial health and economic losses (2). Thus, the study of preventing an epidemic or eradicating an existing epidemic has been a hot topic in the area of public health.

Vaccination is one of the most effective medical methods to control or even eradicate the spread of a sudden epidemic. Over the course of more than two centuries of development, vaccines have helped eliminate epidemics such as smallpox and malaria, and markedly reduced the number of deaths from measles, tetanus and whooping cough each year. However, in the early stage of a vaccine being available, it is impossible to produce sufficient vaccine doses to vaccinate an entire population due to limited vaccine types and production capacity (2). Also, populations with different characteristics, such as age, sex, ethnic group, and occupation, may have different risks of being infectious or dying from a given disease (3, 4); for example, the elderly may have a high probability of becoming infected because most have weakened immune systems. Therefore, allocating limited vaccines over different populations to help more people heal from or avoid the disease entirely is a critical issue in epidemic control.

Many researchers have investigated vaccine allocation since the beginning of the 21st century. Particularly after the COVID-19 outbreak, vaccine development was put on the agenda early in various countries, and major breakthroughs were made quickly (5, 6). Many institutions and research groups began to focus on the optimization of the allocation of the COVID-19 vaccine. These studies can be divided into two types: qualitative and quantitative.

Among the qualitative studies, researchers suggested vaccination priorities based on ethical principles, economic principles, health principles, and so on. In September 2020, the WHO (7) proposed a framework for the allocation and prioritization of COVID-19 vaccination, whose global goal was to let vaccines equitably protect and promote human wellbeing among all people of the world. This framework primarily considered 6 principles: human will-being, equal respect, global equity, national equity, reciprocity and legitimacy. Persad et al. (8) also proposed suggestions, but the ethical values that they emphasized benefitted people and limited harm by prioritizing disadvantaged populations and equally concerning different populations of different genders. Roope et al. (9) focused on socioeconomic development and human health, and raised three different criteria: individual health benefits, social health benefits and economic benefits. A report by the US National Academics of Sciences, Engineering and Medicine (2) summarized the experience of some qualitative research on COVID-19 vaccine distribution and divided the principles to be followed into two categories: ethical principles, and procedural principles. Additionally, the report proposed a four-stage vaccine allocation strategy according to these principles.

Researchers who performed quantitative studies primarily used mathematical models to evaluate the outcomes (e.g., infection and death prevention) of different vaccine allocation strategies and then provided suggestions. Some researchers aimed to optimize vaccine allocation over an entire strategy space (2, 10–14), and others aimed to evaluate and compare pre-given strategies and choose the best one (15–17). Many researchers considered the optimal vaccine allocation over different age groups and thus used an age-stratified SEIR model as an epidemic simulation mathematical model; some also provided deeper considerations. For example, Matrajt et al. (14) and Mylius et al. (17) accounted for the risk of being infected and participated in the population according to age and risk in their models. Both also considered how the allocation time point during an epidemic affected the optimal allocation strategy. Preciado et al. (12) constructed a model at the individual level, used an arbitrary network aggregating an SIS epidemic model to simulate disease transmission in the population, and proposed a convex framework for optimal vaccine allocation over the network. In addition to age, optimal allocation over geographic regions was also studied. Yarmand et al. (13) considered a decision process in which vaccination is performed in two phases and formulated the vaccine allocation problem as a two-stage stochastic linear programming. The optimization objectives or comparison criteria that these studies considered were total infectious, total deaths, years of life lost, economic cost, level of intensive care unit need and so on (11). Even in a recent study of COVID-19 vaccine allocation, a similar research framework was used, for example, the studies of Matrajt et al. (14) and Bubar et al. (16), but the latter assumed that the vaccine may have different effectiveness in the different age groups and thus considered serological tests in their model.

To the best of our knowledge, although many studies have been conducted to solve the vaccine allocation problem when a vaccine shortage exists, and many evaluation criteria have been considered, few studies have considered the trade-off between different criteria, that is, considered multiple criteria simultaneously. In addition, most have considered prevention vaccine allocation or allocation once, but few have considered multistage vaccine allocation during an epidemic. Vaccine allocation is a long, multistage process, and strategies can vary at different stages when vaccine quantities are insufficient; thus, considering multiple stages is critical. The main focus of this research will be developing a framework for determining optimal multistage vaccine allocation strategies with multiple objectives.

In this study, we modeled the vaccine allocation problem as a two-objective optimization problem by considering the total number of disease onsets and the total number of severe cases simultaneously, and solved it with the SPEA2 algorithm. To evaluate the optimization goals, we developed an age-stratified deterministic extended SEIR model to simulate the spread of the disease. We extended the conventional SEIR model to a more complex type by considering more compartments. Using this model, the spread of an infectious disease can be simulated by adjusting certain model parameters. We also considered age-specific contact information and vaccine effects in the extended SEIR model. Additionally, we used a case study to describe the effectiveness of the proposed framework, optimizing the COVID-19 vaccine allocation strategy. Different production capacities, different vaccine effectiveness levels and different vaccine efficiency against death and severity were considered.

Materials and methods

Model formulation

We considered that there are T stages of vaccine allocation and divided the population into M vaccination groups, where M is not greater than the number of age groups (i.e., N). We used X = (x₁₁, x₁₂, …, x_1M, …, x_T1, …, x_TM) to represent a vaccine allocation strategy, whose element x_ti is a representation of the proportion of the number of vaccines allocated to vaccination group i at the t-th vaccine allocation stage to the total number of available vaccines at the allocation stage. Specifically speaking, if the total available vaccine number is M_t, then the number of vaccines allocated to vaccination group i is M_ti = M_t × x_ti. M_t is related to the vaccine production capacity, which is calculated by the product of the vaccine production capacity percent and the population size in our case study. The vaccine allocation strategy among different age groups in the same vaccination group is that the number of vaccines allocated to each age group is proportional to its demands. The two objective functions, f(X) represents the total number of disease onsets, and h(X) represents the total number of severe cases (deaths and cases that require intensive care) under vaccine allocation strategy X. These functions do not have exact expression, and can be evaluated by the age-stratified deterministic extended SEIR model developed in this study. f(X) is the summation of the daily simulated new symptomatic infected individuals; and h(X) is the summation of the sum of the daily simulated new infected individuals requiring ICUs and the daily simulated new deaths. And detailed information about the model and how the strategy X is contained in the vaccination process simulation are described in the section Extended SEIR model. Using these notations, we developed the following optimizing programming model. The constraints mean that the total number of vaccine doses allocated cannot exceed the total number of available vaccines doses at the allocation stage. We used SPEA2 to solve this two-objective optimization problem. The process of SPEA2 and how we use it in this study is described in the section Optimization method:

$\begin{array}{r}{min}_{x}f(\text{X})\\ {min}_{x}h(\text{X})\\ s.t.{\displaystyle \sum _{i=1}^M}{x}_{ti}\le 1,\forall t\\ x_{ti}\ge 0,\forall i,t\end{array}$

Extended SEIR model

Model structure

We developed an age-stratified deterministic extended SEIR model that divides an entire population into N age groups by age. For each age group, as shown in Figure 1, the population is divided into 9 compartments: susceptible (S) population, exposed(E) population, pre-symptomatic (P) population, infectious population, recovered (R) population, and the people who died (D) from the epidemic. Pre-symptomatic individuals can become asymptomatic (A) individuals or symptomatic (I) individuals based on the severity of their symptoms. Those who are symptomatic may recover or die, require hospitalization (H), or require intensive care (C). A capital letter with the subscript V represents the corresponding population who has been vaccinated. The red dotted arrow represents the process of vaccination.

FIGURE 1

Figure 1. Illustration diagram of the extended age-stratified SEIR model.

In this study, due to the short duration of the simulation, we did not consider any births or deaths caused by other reasons and assumed that the vaccines' effect last longer than the time horizon we study. Additionally, we assumed that different age groups had different infection fatality ratios, proportions of hospitalizations requiring ICU admission and susceptibilities. Other characteristics are held constant. Interventions other than vaccines, such as wearing masks and maintaining social distance, were also ignored in this study. We also assumed that one could not be confirmed as an infectious individual except based on for the onset of symptoms. Also, we assumed that at the early stage of vaccination, some testing methods, such as nucleic acid testing, were already mature and freely accessible to all people very easily. This is the real case in China for COVID-19; now almost all residents perform such tests regularly, which gives no chance for asymptomatic to get vaccine. Therefore, considering the scenario in China, we assumed that some tests were performed before one being vaccinated in our model, which can identify the asymptomatic infected individuals. Thus, they have no opportunity to get vaccine. Thus, the populations represented by S, E, and P have the opportunity to get vaccine, and we assumed that one could obtain immunity if and only if he or she gets two vaccine doses.

We described the social contact pattern between the age groups using a social contact matrix ${M}$. Its i-th row and j-th column element m_ij is the expectation value of the number of individuals in age group j that an individual in age group i contacts in a single day. All relevant parameters and their descriptions are shown in Table 1.

TABLE 1

Table 1. Parameter descriptions of the extended age-stratified SEIR model.

Based on these assumptions and descriptions, the model equations without considering vaccination are:

$\begin{array}{l}\frac{d{S}_i}{dt}=-\beta m_i{S}_i{\displaystyle \sum _j}\left(\frac{m_{ij}}{N_j}(r_P(P_j+P_{V,j})+r_A(A_j+A_{V,j})+(I_j+I_{V,j}))\right)\end{array}$

$\begin{array}{l}\frac{d{E}_i}{dt}=\beta m_i{S}_i{\displaystyle \sum _j}\left(\frac{m_{ij}}{N_j}(r_P(P_j+P_{V,j})+r_A(A_j+A_{V,j})+(I_j+I_{V,j}))\right)-\frac{E_i}{D_E}\end{array}$

$\begin{array}{l}\frac{d{P}_i}{dt}=\frac{E_i}{D_E}-\frac{P_i}{D_P}\end{array}$

$\begin{array}{l}\frac{d{A}_i}{dt}=\frac{(1-\alpha )P_i}{D_P}-\frac{A_i}{D_A}\end{array}$

$\begin{array}{l}\frac{d{I}_i}{dt}=\frac{\alpha P_i}{D_P}-\frac{I_i}{D_R}-\frac{I_i}{D_Q}\end{array}$

$\begin{array}{l}\frac{d{C}_i}{dt}=\frac{{\sigma }_i{I}_i}{D_Q}-\frac{C_i}{D_C}\end{array}$

$\begin{array}{l}\frac{d{H}_i}{dt}=\frac{(1-{\sigma }_i)I_i}{D_Q}-\frac{H_i}{D_H}\end{array}$

$\begin{array}{l}\frac{d{R}_i}{dt}=\frac{(1-f_{Hi})I_i}{D_R}+\frac{(1-f_{Ci})C_i}{D_C}+\frac{(1-f_{Hi})H_i}{D_H}+\frac{A_i}{D_A}\end{array}$

$\begin{array}{l}\frac{d{D}_i}{dt}=\frac{f_{Hi}{I}_i}{D_R}+\frac{f_{Ci}{C}_i}{D_C}+\frac{f_{Hi}{H}_i}{D_H}\end{array}$

$\begin{array}{l}i=1,2,\dots ,N\end{array}$

We considered a leaky vaccine that can reduce susceptibility to infection and model it as a reduced probability of getting infectious when a susceptible individual contacts an infectious individual. This capability of reduced susceptibility due to vaccine is referred to as the effectiveness of a vaccine (V_e) in this article. We assumed that V_e were age-independent. According to this assumption, the model equations for describing the population that obtain vaccinated have the same structure as model equations without considering vaccination except for the first and the second equations, and for the vaccinated population, the two different equations are:

$\begin{array}{l}\frac{d{S}_{V,i}}{dt}=-\beta (1-V_e)m_i{S}_{V,i}{\displaystyle \sum _j}\left(\frac{m_{ij}}{N_j}(r_P(P_j+P_{V,j})+r_A(A_j+A_{V,j})+(I_j+I_{V,j}))\right)\end{array}$

$\begin{array}{l}\frac{d{E}_{V,i}}{dt}=\beta (1-V_e)m_i{S}_{V,i}{\displaystyle \sum _j}\left(\frac{m_{ij}}{N_j}(r_P(P_j+P_{V,j})+r_A(A_j+A_{V,j})+(I_j+I_{V,j}))\right)-\frac{E_{V,i}}{D_E}\end{array}$

For the process of vaccination, we assumed that a susceptible individual acquired immunity once he or she is vaccinated, and then we can use the following equations to model the process:

$\begin{array}{r}{S}_{V,i,t_1}=S_{V,i,{(t-1)}_L}+min\left\{S_{i,{(t-1)}_L},\frac{S_{i{,(t-1)}_L}}{S_{{i,(t-1)}_L}+E_{i,{(t-1)}_L}+P_{{i,(t-1)}_L}}{M}_{ti}\right\}\\ E_{V,i,t_1}=E_{V,i,{(t-1)}_L}+min\left\{E_{i,{(t-1)}_L},\frac{E_{i{,(t-1)}_L}}{S_{{i,(t-1)}_L}+E_{i,{(t-1)}_L}+P_{{i,(t-1)}_L}}{M}_{ti}\right\}\\ P_{V,i,t_1}=P_{V,i,{(t-1)}_L}+min\left\{P_{i,{(t-1)}_L},\frac{P_{i{,(t-1)}_L}}{S_{{i,(t-1)}_L}+E_{i,{(t-1)}_L}+P_{{i,(t-1)}_L}}{M}_{ti}\right\}\\ S_{i,t_1}=S_{i,{(t-1)}_L}-min\left\{S_{i,{(t-1)}_L},\frac{S_{i{,(t-1)}_L}}{S_{{i,(t-1)}_L}+E_{i,{(t-1)}_L}+P_{{i,(t-1)}_L}}{M}_{ti}\right\}\\ E_{i,t_1}=E_{i,{(t-1)}_L}-min\left\{E_{i,{(t-1)}_L},\frac{E_{i{,(t-1)}_L}}{S_{{i,(t-1)}_L}+E_{i,{(t-1)}_L}+P_{{i,(t-1)}_L}}{M}_{ti}\right\}\\ P_{i,t_1}=P_{i,{(t-1)}_L}-min\left\{P_{i,{(t-1)}_L},\frac{P_{i{,(t-1)}_L}}{S_{{i,(t-1)}_L}+E_{i,{(t-1)}_L}+P_{{i,(t-1)}_L}}{M}_{ti}\right\}\end{array}$

where M_ti is the number of people who can obtain vaccinated using the vaccines allocated to age group i in the t-th time vaccine allocation, which is related to the vaccine production capacity, with two vaccine doses per person; L is the time lag between two adjacent vaccine allocation stages; t₁ is the time point after the t-th time vaccine allocation; and (t − 1)_L is the time point before the t-th time vaccine allocation.

We take the first and the fourth equations as examples to explain the meaning of these eight equations. The first equation represents the immediate change of the size of vaccinated susceptible population i after the t-th time vaccine allocation. S_V,i,(t−1)L is the size of vaccinated susceptible population i at the moment before t-th time vaccine allocation, and S_{V, i,t1} is the size of vaccinated susceptible population i immediately after t-th time vaccine allocation is completed. Under the assumption that the number of vaccines allocated to one compartment group is proportional to its size, $min\left\{S_{i,{(t-1)}_L},\frac{S_{i{,(t-1)}_L}}{S_{{i,(t-1)}_L}+E_{i,{(t-1)}_L}+P_{{i,(t-1)}_L}+A_{{i,(t-1)}_L}}{M}_{ti}\right\}$ is the number of vaccines allocated to susceptible population i in the t-th time vaccine allocation. min{a, b} represents the smaller one between a and b, we used this notation to prevent the number of vaccines allocated to one compartment group to exceed its needs. And the fourth equation represents the immediate change of the size of unvaccinated susceptible population i after the t-th time vaccine allocation. Similarly, the other six equations represent the immediate change of the other compartment groups after the t-th time vaccine allocation.

Transmission coefficient calculation

Most disease parameters in the model, such as the mean duration of the latent period, can be determined from existing research about the disease, but the transmission coefficient β is not available. Due to policies and associated regularities made by the government, the social activity patterns of the population tend to change during an epidemic, which leads to a variational social contact matrix ${M}$; thus, we cannot model this variational pattern due to its complexity, and using real epidemic data to estimate β is nearly impossible.

This study aims to optimize the multi-time vaccine allocation strategy without considering other interventions; thus, we assumed that the ${M}$ during an epidemic remained the same as it is in the normal situation, and the ${M}$ in the normal situation was available because a few studies have given their estimation of ${M}$ using some statistical methods based on some investigation data (18). We can thus use their results in this study. Many researchers have estimated the basic reproduction number R₀ in the early stage of the epidemic due to its critical role in infectious disease studies, such as the R₀ of COVID-19 (22). R₀ has a tight relationship with nearly all the model parameters, including ${M}$ and β. This reality inspired us to use all the other parameter values to calculate β directly. The next generation matrix (NGM) provided us with such a possibility. Using NGM to determine R₀ for an ordinary differential equations (ODE) model was proposed by Diekmann et al. (23) and elaborated by Van den Driessche and Watmough (24). Thus, NGM helps us to construct a mathematical relation between R₀ and model parameters. The following provides a brief introduction about how to construct an NGM and how to use it to calculate R₀ and further displays the method of β calculation.

We assumed that the entire population was separated into n compartments and let x = (x₁, x₂, …, x_n) be the number of individuals in each compartment. Without a loss of generality, we assumed that the first m < n compartments contained infectious individuals; then, the ODE of the first m compartments can be rewritten in the form of $\frac{d{x}_i}{dt}={{F}}_i(\text{x})-{{V}}_i(\text{x})$, where ${{F}}_i(\text{x})$ is the rate of appearance of new infections in compartment i, and ${{V}}_i(\text{x})$ is the rate of other transitions between compartment i and other infected compartments. Using the concept of the Jacobian matrix, we defined $F={\left[\frac{d{F}_i({\text{x}}_0)}{d{x}_j}\right]}_{m\times m},V={\left[\frac{d{V}_i({\text{x}}_0)}{d{x}_j}\right]}_{m\times m}$, where x₀ is a disease freedom equilibrium (DFE) point of the system. Then, FV⁻¹ is the NGM of the model, and R₀ is the spectral radius of the NGM.

In the proposed model, we assumed that compartments E, P, A, and I contained infectious individuals. According to the splitting principles mentioned above, the ${F}$ functions, ${V}$ functions, and DFE should be:

$\begin{array}{l}{F}=\left(\frac{\beta m_i{S}_i{\displaystyle {\sum }_j}\frac{m_{ij}}{N_j}(r_P{P}_j+r_A{A}_j+I_j)}{\begin{array}{c}0\\ 0\\ 0\end{array}}\right),i=1,2,\dots ,N\end{array}$

$\begin{array}{l}{V}=\left(\begin{array}{c}\frac{E_i}{D_E}\\ \frac{P_i}{D_P}-\frac{E_i}{D_E}\\ \frac{A_i}{D_A}-\frac{(1-\alpha )P_i}{D_P}\\ \frac{I_i}{D_R}+\frac{I_i}{D_Q}-\frac{\alpha P_i}{D_P}\end{array}\right),i=1,2,\dots ,N\end{array}$

$\begin{array}{l}{x}_0=(S_{0,1\times N},E_{0,1\times N},P_{0,1\times N},A_{0,1\times N},I_{0,1\times N},H_{0,1\times N},)\\ {(C_{0,1\times N},R_{0,1\times N},D_{0,1\times N})}^T\\ ={(N_1,N_2,\dots ,N_N,0_{1\times 8N})}^T\end{array}$

Defining ${{M}}^{{}^{\prime }}={\left[m_{ij}^{{}^{\prime }}\right]}_{N\times N}={\left[\frac{m_i{m}_{ij}{S}_{i0}}{N_j}\right]}_{N\times N}$, we provide the formulation of F and V as follows:

$\begin{array}{l}F=\left[\begin{array}{c}{0}_{N\times N}\beta r_P{{M}}^{{}^{\prime }}\beta r_A{{M}}^{{}^{\prime }}{\beta {M}}^{{}^{\prime }}\\ 0_{N\times N}{0}_{N\times N}{0}_{N\times N}{0}_{N\times N}\\ 0_{N\times N}{0}_{N\times N}{0}_{N\times N}{0}_{N\times N}\\ 0_{N\times N}{0}_{N\times N}{0}_{N\times N}{0}_{N\times N}\end{array}\right]\\ =\beta \left[\begin{array}{c}{0}_{N\times N}{r}_P{{M}}^{{}^{\prime }}{r}_A{{M}}^{{}^{\prime }}{{M}}^{{}^{\prime }}\\ 0_{N\times N}{0}_{N\times N}{0}_{N\times N}{0}_{N\times N}\\ 0_{N\times N}{0}_{N\times N}{0}_{N\times N}{0}_{N\times N}\\ 0_{N\times N}{0}_{N\times N}{0}_{N\times N}{0}_{N\times N}\end{array}\right]=\beta F^{{}^{\prime }}\end{array}$

$\begin{array}{l}V=\left[\begin{array}{c}\frac{1}{D_E}{I}_{N\times N}{0}_{N\times N}{0}_{N\times N}{0}_{N\times N}\\ -\frac{1}{D_E}{I}_{N\times N}\frac{1}{D_P}{I}_{N\times N}{0}_{N\times N}{0}_{N\times N}\\ 0_{N\times N}-\frac{(1-\alpha )}{D_P}{I}_{N\times N}\frac{1}{D_A}{I}_{N\times N}{0}_{N\times N}\\ 0_{N\times N}-\frac{\alpha }{D_P}{I}_{N\times N}{0}_{N\times N}\left(\frac{1}{D_R}+\frac{1}{D_Q}\right)I_{N\times N}\end{array}\right]\end{array}$

where 0_{N × N} is a zero matrix with dimensions of N × N, and I_{N × N} is an identity matrix with dimensions of N × N. We used ρ to represent the spectral radius of a matrix. Then, the following mathematical relation yields us a method to calculate β by other model parameters and R₀:

$\begin{array}{l}{R}_0=\rho \left(F{V}^{-1}\right)=\beta \rho \left(F^{{}^{\prime }}{V}^{-1}\right)\end{array}$

Model extension

We discuss an important problem in this subsection, which may significantly influence the optimal vaccine allocation results. Some vaccines not just protect one individual from being infected, but reduce the infected fatality rate and the proportion of the hospitalization individuals requiring ICU. To tackle this important issue, we further extended our SEIR model, the model structure is shown in Figure 2. To represent the vaccine efficiency against death and severity, we added a parameter V_ep in the SEIR model. Then the infected fatality rate and the proportion of hospitalization individuals requiring ICU of vaccinated group are $f_{C(H)i}^V={(1-V_{ep})f}_{C(H)i}$ and ${\sigma }_i^V=(1-V_{ep}){\sigma }_i$, respectively. The model equations for describing the population that obtain vaccinated in this scenario have the same structure as model equations considering vaccination in section Model structure except for replacing f_C(H)i and σ_i with $f_{C(H)i}^V$ and ${\sigma }_i^V$.

FIGURE 2

Figure 2. Illustration diagram of the further extended age-stratified SEIR model.

Optimization method

In this study, we considered two optimization objectives simultaneously: the total number of onsets and the total number of severe results (cases that require intensive care and deaths). We then optimized the vaccine allocation strategy with different vaccine effectivenesses, production capacities and efficiency against death and severity.

We must solve a two-objective optimization problem. The methods used to solve a multi-objective problem can be divided into two types (25, 26). The first type is to transform the original optimization problem into a single objective optimization problem and then to solve the new problem by some single-objective optimization method; however, these methods may have some strong restrictions on the property of objective space. The second type is to solve the original problem directly and to find a set of optimal solutions named Pareto optimal solutions; however, these methods cannot typically find all true Pareto optimal solutions and only obtain a set of nondomination solutions as an approximation. Almost all of these methods are evolutionary methods and are used widely in the area of multi-objective optimization; however, these methods are time-consuming. In this study, we used the Improved the Strength Pareto Evolutionary Algorithm (SPEA2), which was proposed by Zitzler et al. (27). Based on the Strength Pareto Evolutionary Algorithm (SPEA), Algorithm 1 shows the framework of SPEA2, and more details of SPEA2 are available in the reference paper.

ALGORITHM 1

Algorithm 1. SPEA2 framework.

An individual of SPEA2 has the same structure as the allocation strategy X (decision variables) that we defined in the section Model formulation, and each element of the vector X is a gene of the related individual; thus, the individual is real-coded in this problem. Thus, we used extended intermediate crossover (28) and Mühlenbein's mutation (28) in Step 6.

Note that SPEA2 cannot guarantee the feasibility of a solution, and the programming problem that we formulated in the Section Model formulation has two types of constraints. Infeasible solutions may be generated in Step 6 that must be ‘repaired’. First, if there are some x_ijs in an individual being < 0, we must generate a random number from a uniform distribution between 0 and 1 to replace it. Then, if a solution X with $\sum _{i=1}^N{x}_{ti}>1$ for some t, we must map it to X' such that $\sum _i^N{x}_{ti}^{{}^{\prime }}>1$ for those t by replacing x_ti with $x_{ti}^{{}^{\prime }}=x_{ti}/\sum _{i=1}^N{x}_{ti}$.

There may also be another type of potentially infeasible situation. Even if a solution meets the constraints of the programming model, it may mislead the analysis and may yield confusing or incorrect suggestions. This type of situation is where vaccines that are allocated to some age groups at some vaccine allocation times exceed their total demands at that time; thus, a portion of vaccines are wasted. By replacing the proportion with the true demand proportion, we can ‘repair’ such an infeasible situation.

Results

Simulation settings

We applied the framework proposed in this study to COVID-19 vaccine allocation in a population whose age structure is similar to Wuhan City, China. The COVID-19 epidemic has led to tens of millions of deaths and heavily impeded economic development (29). A number of modelers have made efforts to develop a variety of models to study the spread of the epidemic (19, 30–39). These studies have provided adequate information on the characteristics of the infectious disease, such as the basic reproduction number and duration of the latent period (22). Most parameter values of the simulation model we used in this section were borrowed from previous studies (Table 1). It should be noted that the infected fatality rate of individuals requiring ICUs is greater than that of not requiring ICUs. However, we treated deaths and cases requiring ICUs both as severe cases, thus, the death rate for individuals requiring ICUs is not so important in this study. In addition, the overall death rate approximately equals the death rate of individuals requiring hospitalization due to the relatively low severity rate of COVID-19 (1). Hence, we assumed f_Ci = f_Hi = f_i, and used the overall death rate to approximate the death rate of individuals requiring hospitalization in our case study, which is reasonable.

Since we wanted to simulate the spread of COVID-19 among different age groups without any other interventions except vaccine injection, a social contact matrix ${M}$ was required, which can describe the contact pattern of the population mentioned above in the normal situation. We used the ${M}$ obtained by Zhang et al. (18) because it described the social contact patterns of citizens in Wuhan from December 24, 2019, to December 30, 2019, when citizens of Wuhan City had no sense of the outbreak of the epidemic. ${M}$ divides the entire population into 14 age groups (N = 14): 0–4, 5–9, 10–14, 15–19, 20–24, 25–29, 30–34, 35–39, 40–44, 45–49, 50–54, 55–59, 60–64, and 65+; Figure 3 shows a heatmap of ${M}$. For the 65+ age group, we performed a weighted average of the rates in the Table 2 according to the relative percentages of the population aged 65–69, 70–74, 75–79, and 80+.

FIGURE 3

Figure 3. Heatmap of ${M}$ in Wuhan City.

TABLE 2

Table 2. Values of N_i, f_i, σ_i, and m_i of different age groups.

We obtained the population size of each age group in Hubei Province from the Sixth Nation Population Census in 2010 from the official website of the National Bureau of Statistics of China (40) and estimated the age structure of Wuhan citizens from the data. Aggregating with the population size of Wuhan in 2020, we calculated the size of each age group in the population that we considered to perform simulations while ignoring the age structure changes during the last decade. Note that the basic reproduction number we used to calculate β was estimated using the data from the early stage of the epidemic. We set D_Q = 21 days when calculating β, and we set D_Q = 1.5 days when performing simulation because the epidemic prevention and control is strict now, and the time lag between symptom onset and hospitalization is shorter than 2 days.

Figure 3 shows that in the normal situation, the highest social contacts occur between 5 and 19 age groups and their peers. The contact between 20 and 49 age group and their peers is also frequent, but the frequency is lower than that between 5 and 19 age groups and their peers. Individuals in 20–49 age group also frequently contact the 50–64 age group. The 50–64 age group contacts both the 20–49 age group and the 65+ age group frequently, which indicates that 20–49 and 50–64 age groups may have a higher transmission risk. This analysis inspired us to consider 4 vaccination groups (M = 4): 0–19 age group, 20–49 age group, 50–64 age group and 65+ age group. We considered 3 vaccine allocation stages (T = 3), and the time lag between any two adjacent allocation time points was 30 days (L = 30). To investigate how vaccine effectiveness and vaccine production capacity affect optimization, we considered 9 vaccine effectiveness (V_e from 10 to 90% in steps of 10%) and 5 vaccine production capacities (VPCs) (5, 10, 15, 20, 25%). Vaccine production capacity was defined as the proportion of the number of vaccine doses available at each allocation time point over the size of the entire population that we considered, and we assumed that the capacity was the same at different stages. We only set V_ep = 0, 0.3, 0.6, 0.9 when V_e = 0.3, 0.6, 0.9 in the model extension. In this section, we assumed that all individuals were susceptible; thus, the size of the population with pre-existing immunity was zero, and in real applications, this setting could be changed when different real scenarios are met.

The situation we considered in this section was an imported epidemic, where several infected individuals are introduced into an epidemic-free area and lead to an epidemic outbreak. Thus, the initial simulation state we set was one infected individual in each age group, and the others were all susceptible individuals. The time length we considered to calculate the total onsets and severe cases was 200 days. Additionally, the baseline strategy we used to compare with optimal strategies was that the number of vaccines allocated to each age group is proportional to its demands. Under the above settings, the optimal allocation strategies were obtained using SPEA2 algorithm descripted in section Optimization method.

Optimization results analysis

To display and analysis the final optimization results, we introduce some basic concepts in the field of multi-objective optimization (41):

• Pareto domination: A solution X dominates another solution Y if all the objectives of X are not worse than Y's, and X is strictly better than Y in at least one objective.

• Pareto optimality: A solution X is a Pareto optimal solution if there is no solution in the feasible region dominating X.

• Global Pareto optimal set: A global Pareto optimal set consists of all the Pareto optimal solutions in the feasible region.

• Pareto optimal front: The projections of all the elements of the global Pareto optimal set in the objective space constitute the Pareto optimal front.

Solving a multi-objective problem requires identifying its Pareto optimal front; however, finding the true Pareto optimal front is generally impossible. Thus, the primary aim of solving a multi-objective optimization problem is to find a well-distributed, approximate Pareto optimal front. A Pareto figure is a tool to display the optimization results in the field of multi-objective optimization. Usually, the two axises represent different objective function values; points on the figure represents the projections of feasible (or Pareto optimal) solutions in the objective space.

We drew Pareto figures under the different combinations of vaccine effectiveness and vaccine production capacity (Figure 4; Supplementary Figures S1–S8). The subfigures A–E in Figures 4, 5 are the Pareto figures under different vaccine production capacities, and the relevant vaccine production capacity is displayed below each subfigure. In each subfigure, the abscissa axis is the total severe cases (h(X)), and the ordinate axis is the total onsets (f(X)). Different subfigures in each figure show the results under different vaccine production capacities. The black points in each figure are the simulation results of 800 random strategies; the blue points are the approximate Pareto optimal front obtained by SPEA2, the simulation outputs of the Pareto optimal solutions; and the red point is the simulation output of the baseline strategy. We could easily notice that the black points and the red points always lie on the upper right side of the approximate Pareto optimal front, which means that the optimization process is important. Also, in the scope of this study, along with the increase in vaccine effectiveness and vaccine production capacity, we discovered that the contradiction between the two objectives becomes weaker, and the positive correlation property becomes stronger. Consequently, the Pareto optimal front is broad in a low vaccine effectiveness (not > 30%, Figure 4; Supplementary Figures S1, S2) situation and narrow with relatively high vaccine effectiveness (not < 70%, Supplementary Figures S6–S8), while with a medium vaccine effectiveness (say, between 30 and 70%, Supplementary Figures S3–S5), the Pareto optimal front is broad under a low vaccine production capacity and narrow under a relatively high vaccine production capacity.

FIGURE 4

Figure 4. Pareto optimal fronts when V_e = 10%.

FIGURE 5

Figure 5. Heatmaps of Pareto optimal solutions when V_e = 10%.

In an effort to explore the detailed structure of the Pareto solutions and how the vaccine parameters affect the Pareto solutions, we drew a heatmap of the Pareto solutions under each combination of vaccine parameters (Figure 5; Supplementary Figures S9–S16). The subfigures A–E in Figures 4, 5 are the heatmap of the Pareto solutions under different vaccine production capacities, and the relevant vaccine production capacity is displayed below each subfigure. Each row of each heatmap is a corresponding Pareto solution, and the number of its left side is the solution id. The solution numbered ‘baseline’ is the baseline solution. In each heatmap, the direction of solution id increasing is the direction of the total onsets increasing but total severe cases decreasing. The color of each square is the related value of x_ti.

The structure of the baseline solution indicates that the baseline strategy always allocates most vaccines to the 20–49 age group, allocates the second most vaccines to the 0–19 age group, and allocates the fewest number of vaccines to the 65+ age group, regardless of vaccine parameters and allocation stage when the vaccine is insufficient. This baseline strategy strongly positively correlates with the age group size, which is consistent with our perceptual intuition.

Based on the heatmaps of the Pareto optimal solutions, we summarized the rules described in the optimal allocation strategy under different vaccine effectiveness and vaccine production capacity. These rules were only discussed within the value range of each parameter considered in this case study.

If we look at the 0–49 age group as a whole, and the 50+ age group as a whole, and consider the three vaccine allocation phases together, the optimal strategies that tend to minimize the total number of onsets suggest allocating most of the vaccines to the 0–49 age group, while optimal strategies that tend to minimize the total number of severe results (i.e., deaths or ICU cases) suggest allocating most of the vaccines to the 50+ group. This tendency diminishes as vaccine effectiveness increases and vaccine production capacity increases. When the vaccine effectiveness and the vaccine production capacity reach a certain level, the contradiction between the two goals is weakened, the positive correlation is enhanced, the Pareto front becomes narrow, we can only obtain a small number of Pareto optimal solutions, and this trend tends to disappear.

Note that people aged 65+ have a higher risk of being infected and dying from the infectious, and we referred the strategy that directly allocates most of the vaccines to this population as the direct protection of a high-risk population. Conversely, although young and middle-aged groups have a lower rate of infection and death, they frequently contact most age groups and thus have a higher transmission risk. Protecting these age groups can control epidemic spread, and we thus referred the strategies that directly protect these age groups as the indirect protection of high-risk populations.

Based on this analysis, we provided a more detailed discussion. When the vaccine effectiveness is low (10, 20, and 30%, Figure 5; Supplementary Figures S9, S10), and with increasing vaccine production capacity, the vaccines allocated to the 0–49 age group gradually shift from being more allocated to the 0–19 age group to being more allocated to the 20–49 age group. Also, the vaccines allocated to the 50+ age group gradually shift from being more allocated to the 65+ age group to the 50–64 age group (i.e., shifting from direct protection to indirect protection). The results in each allocation stage are consistent. When the vaccine effectiveness is relatively high (60, 70, 80, and 90%, Supplementary Figures S13–S16), and the vaccine production capacity is low (5%), the optimal strategy suggests allocating most of the vaccine to the 65+ age group to directly protect them. As vaccine production capacity increases, the optimal strategy suggests allocating most of the vaccines to the 20–49 and 50–64 age groups. The strategy with a higher vaccine effectiveness in this range tends to allocate more vaccines to the 20–49 age group and fewer to the 50–64 age group. This trend is enhanced when the allocation stage increases. Additionally, with medium vaccine effectiveness (40 and 50%, Supplementary Figures S11, S12), and when vaccine production capacity is low (5%), the optimal strategy allocates most vaccines to the 20–49 or 65+ age group according to the optimizing tendency, and the results at each allocation stage are consistent. As vaccine production capacity increases, the optimal strategy allocates a proportion of vaccines to the 65+ age group, which decreases as the allocation stage increases and then primarily allocates the vaccines to the 20–49 or 50–64 age group in line with the optimization tendency. When the vaccine effectiveness and vaccine production capacity are sufficiently high (90%, 25%, Supplementary Figure S16E), the optimal strategy primarily allocates vaccines to 20–49 and 65+ age groups, and with increasing allocation, the proportion allocated to 20–49 age group increase, while the proportion allocated to 65+ age group decrease.

Simulation comparation of optimal strategy and baseline strategy

To evaluate the optimal strategy in more detail, we chose one solution for each combination of vaccine parameters to perform a simulation analysis and compare the simulation outcomes with the outcomes of the baseline strategy. Note that the difference in the magnitude of the two objectives is large. We scaled them to values between 0 and 1 using Equations (2) and (3), where D is the nondomination set obtained under a combination of vaccine parameters, and then choose the solution that is the closest to the original point after scaling (described as Equation 1) in the objective space. As shown in the Pareto figures, with some vaccine parameters, the number of Pareto solutions is not > 2 due to the positive correlation property between the two objectives. In this scenario, we chose one Pareto optimal solution randomly:

$\begin{array}{ll}{\displaystyle \underset{{\text{X}}_i\in \text{D}}{min}}\sqrt{{\left(h^s\left({\text{X}}_i\right)\right)}^2+{\left(f^s\left({\text{X}}_i\right)\right)}^2}& (1)\end{array}$

$\begin{array}{ll}{h}^s\left({\text{X}}_i\right)=\frac{h\left({\text{X}}_i\right)-\text{ min}\left\{h\left({\text{X}}_j\right)|{\text{X}}_j\in \text{D}\right\}}{max\left\{h\left({\text{X}}_j\right)|{\text{X}}_j\in \text{D}\right\}-\text{ min}\left\{h\left({\text{X}}_j\right)|{\text{X}}_j\in \text{D}\right\}}& (2)\end{array}$

$\begin{array}{ll}{f}^s\left({\text{X}}_i\right)=\frac{f\left({\text{X}}_i\right)-\text{ min}\left\{f\left({\text{X}}_j\right)|{\text{X}}_j\in \text{D}\right\}}{max\left\{f\left({\text{X}}_j\right)|{\text{X}}_j\in \text{D}\right\}-\text{ min}\left\{f\left({\text{X}}_j\right)|{\text{X}}_j\in \text{D}\right\}}& (3)\end{array}$

Figure 6 shows the simulation results of the baseline strategy and the chosen Pareto optimal strategy under different vaccine parameter combinations. Different subfigures are the results under different vaccine effectiveness, and relevant vaccine effectiveness is shown in Figures 6A–I. The abscissa axis is the simulation time, and the ordinate axis is the daily onset number. In each subfigure, the black solid line is the simulation result without a vaccine, the lines in other colors are the simulation results with a vaccine, and different colors stand for different vaccine production capacities. Solid lines show the baseline strategy simulation results, and dashed lines show the chosen Pareto optimal strategy simulation results. There are two primary findings from these simulations. First, vaccination can delay the arrival of the peak daily onset number and suppress it. Additionally, the ability of vaccination to delay and suppress the peak is enhanced as the two vaccine parameters increase. Second, the ability of the Pareto optimal strategy to delay and suppress the peak is more permanent than that of the baseline strategy.

FIGURE 6

Figure 6. Simulation comparation between the chosen optimal strategy and baseline.

Also, considering the uncertainty of some model parameters, we randomized some model parameters (R₀, α, D_E, D_P, D_C, D_H, D_R, and D_Q) and assigned a continuous probability distribution to each of these parameters. Then, we sampled 1,000 groups of parameters from these distributions. Finally, we ran 1,000 simulations using these parameters. Figure 7 and Supplementary Figures S17, S18 compare the chosen Pareto optimal strategy and baseline strategy with regard to onset prevention, ICU case prevention and death prevention separately under different combinations of vaccine parameters, respectively. Different subfigures show the results under different vaccine effectiveness, and relevant vaccine effectiveness is in Figures 6A–I. The abscissa axis is the vaccine production capacity, and the ordinate axis is the percentage of onset prevention, ICU case prevention or death prevention. The blue line shows the results of the chosen Pareto optimal strategy, and the green line shows the results of the baseline strategy. The shaded area of the corresponding color shows the range of 1,000 simulations without considering the least 2.5% and the most 2.5% simulations.

FIGURE 7

Figure 7. Onsets prevention comparation between chosen optimal strategy and baseline.

The figures indicate that both the baseline strategy and the chosen optimal strategy can prevent onsets, ICU cases and deaths, even if the least effective vaccine (10%) and the least production capacity (5%) are considered. The chosen Pareto optimal strategy is always better than the baseline strategy in all three evaluations. As the two vaccine parameters increase, the effectiveness of vaccines in preventing onsets, ICU cases and deaths becomes increasingly strong, and the advantages of the chosen Pareto optimal strategies become increasingly strong to some extent. When the vaccine is 80% effective and the vaccine production capacity exceeds 20%, the percentage of onset prevention, ICU case prevention and death prevention of the chosen strategy approaches or exceeds 90%. The range between the 2.5% limit and the 97.5% limit of 1,000 simulations is large, which indicates that the effectiveness of the Pareto optimal strategy is sensitive to the model parameters. Although the blue and green shaded areas overlap, the 2.5 and 97.5% limits of the chosen Pareto optimal strategy are nearly always greater than those of the baseline strategy separately, which highlights the advantages of the chosen optimal strategies.

Optimization results analysis of extended model

We accounted for the vaccine that can reduce the infected fatality rate and the proportion that the hospitalization individuals requiring ICU. As is shown in section Model extension, we used a parameter V_ep to model this effect. We reoptimized the allocation policies under different combinations of V_e (30, 60, 90%), vaccine production capacity (VPC) (5, 10, 15, 20, 25%), and V_ep (0.3, 0.6, 0.9), analyzed the results and compared them with the results obtained when V_ep = 0 (i.e., did not consider the influence of vaccine on the infected fatality rate and the proportion that the hospitalization individuals requiring ICU).

Supplementary Figures S19–S21 show Pareto optimal fronts under different combinations of V_e, VPC, and V_ep. Each figure displays the results of a value of V_e; each subfigure displays the results of a value of VPC; and in each subfigure, scatter points of different shapes represent results under different V_ep. For the results under each combination of parameters, the yellow scatter points in each figure are the simulation results of 800 random strategies; the blue scatter points are the approximate Pareto optimal front obtained by SPEA2, the simulation outputs of the Pareto optimal solutions; and the red scatter point is the simulation output of the baseline strategy. We could find that the total onsets are the same but the total severe cases are different for a given vaccine allocation policy, vaccine with a higher V_ep has a lower number of the total severe cases. Consequently, for each combination of VPC and V_ep, the Pareto optimal front with a higher V_ep is lower in the vertical direction than the Pareto front with a lower V_ep, but there is little difference in the horizontal direction. As VPC and V_e increasing, this effect is weakened, and these Pareto fronts get so close to each other that we cannot distinguish them from each other in the picture for high enough VPC and V_e (Supplementary Figure S21E). These findings may be explained by the following two reasons. Firstly, we assumed that the recovered, dead, and hospitalized individuals had no potential to infect others, and V_ep can only affect the size of these three populations. In other words, the infected pressure of susceptible population is the same under different V_eps. Thus, V_ep has no effect on the total onsets, but has an effect on the total severe cases under the same vaccine allocation strategy. Secondly, the optimal policies under a high VPC and V_e suggest to allocate most vaccine doses to the high-transmission population no matter what value of V_ep is (the last row of Supplementary Figure S24), whose infected fatality rate and proportion of the hospitalization requiring ICU are low. Therefore, most of the severe cases are from high-risk population, and changing V_ep have a slight influence on it.

Supplementary Figures S22–S24 show heatmaps of Pareto optimal solutions of model extension. Each figure displays the results of a value of V_e, each row in a figure shows the results of the same VPC and different V_eps, and each column in a figure shows the results of the same V_ep and different VPCs. We could find that the patterns of optimal allocation policies under different V_eps (each row of each figure) are similar, which indicates that V_ep does not have any significant influence on the optimal allocation policies in the scope of this study.

Discussion

This study demonstrated the use of an extended deterministic age-stratified modeling approach together with a multi-objective optimization framework to optimizethe vaccine allocation strategies, and considered the trade-off between the two objectives, i.e., minimizing total onset cases and minimizing total severe cases. We mainly accounted for how different vaccine efficiencies and different vaccine production capacities influence on the optimal policy.

We performed a case study of COVID-19 vaccine allocation to validate the efficiency of the proposed method. One basic finding is that the optimal policies tend to minimize total severe cases prioritize the older population, while the optimal policies tend to minimize total onset cases prioritize the younger population. This finding is in agreement with some previous works that evaluate or optimize the COVID-19 vaccine allocation policies (14, 16, 42–45). This is mainly because the infected fatality rate of COVID-19 is highly positive correlation with age, that is, the older one typically has a higher fatality rate. However, Sunohara et al. (46) suggested that prioritized younger generation was better in terms of deaths under an assumption of a linear relationship between lockdown intensity and acceptable economic loss, and that under an assumption of non-linear relationship, the old first policies were best considering small basic reproduction number. We believe that the reason for the difference in results is that we did not consider factors such as lockdown and economic effect, but only considered the effect of vaccine effectiveness and productivity on the optimization results. To some extent, this shows that in order to obtain more reasonable optimization results, more factors need to be considered according to the actual situation. The results also indicate that the conflict between the two objectives is weakened as the increasement of the vaccine efficiency and the vaccine production capacity because of the increasement of vaccine coverage. What's more, vaccine efficiency and vaccine production capacity are able to influence the switch pattern of priority age group when the optimization tendency shifts from minimizing total onsets to minimizing total severe cases according to the results. For instance, with a low vaccine efficiency and a low vaccine production capacity, the priority age group switches from 0–19 to 65+ when the optimization tendency shifts from minimizing total onsets to minimizing total severe cases; but with a low vaccine efficiency and a high vaccine production capacity, the priority age group switches from 20–49 to 50–64 when optimization tendency shifts. In reality, the other factor, such as the vaccine type (leaky vaccine or all-or-nothing vaccine), can also affect the switch pattern (16). Although we did not consider any other factors except for vaccine efficiency and vaccine production capacity, we provided a detailed description for how these two factors influence priority switch pattern.

The policies we considered in the study are dynamic, that is, the proportion of vaccines allocated to each age group changes as the epidemic progresses. Some researchers also discussed this type of policies in their studies (42, 45). These studies argued that high priority groups—by percent of group vaccinated—switched from one to another as the epidemic develops. Due to the limitations of the optimization method, we only considered three allocation periods, so our results did not show similar phenomena. But the dynamic policies can definitely improve the efficiency of vaccine use due to its flexibility (42).

In summary, within the scope of this study, when vaccine effectiveness is low, the better policy is to allocate most vaccines to the 0–19 age group or 65+ age group under a low production capacity according to the optimization tendency and allocate most vaccines to the 20–49 age group or 50–64 age group according to the optimization tendency under a relatively high production capacity. When vaccine effectiveness is high, the better policy is to allocate vaccines to the 65+ age group with a low production capacity and allocate vaccines to the 20–49 age group with a relatively high production capacity. Matrajt et al. (14) found that under high vaccine coverage, the optimal allocation for all objective functions shifted toward the high-transmission groups. If we treat 20–49 and 50–64 age groups as high-transmission groups, we obtained the similar pattern. However, Matrajt et al. (14) considered one-time vaccine allocation, while we considered multi-period vaccine allocation. In addition, considering the vaccine that can reduce infected fatality rate and proportion of hospitalization requiring ICU do not significantly affect the optimal policy, but can reduce total severe cases.

In the parameter uncertainty study, the 95% confidence intervals are wide, which indicates that the optimal policies may be sensitive to these key parameters. This inspires us to optimize vaccine allocation policy using as precise parameters as possible. Moreover, the existence of other Non-medical interventions may influences the contact pattern, and further influences optimization results, such as the results of Sunohara et al. (46) which considered lockdown.

Conclusion

When an epidemic breaks out, a vaccine, if available, is one of the most useful and effective tools to control and even eliminate the epidemic. With the development of science and technology, the speed of vaccine research and development is becoming increasingly fast. However, in the early stage of vaccine production, the production capacity is low, and the vaccine supply is limited. Thus, policy-makers must decide how to allocate insufficient vaccines among different age groups or risk groups. Many researchers have investigated this important issue; however, most considered different optimization objectives and modeled vaccine allocation optimization as a single objective optimization problem but ignored the trade-off among different optimization objectives.

Acknowledging the shortcomings of past studies, we proposed a new framework to solve the problem of vaccine allocation optimization. In this framework, our contribution is three-fold. Firstly, we consider allocating vaccines among different age groups and modeled this process as a two-objective optimization problem to consider the two optimization objectives simultaneously and discover the trade-off relationship between them. If necessary, the model can be extended to more than two optimization objectives. Secondly, multistage allocation was also considered to model the scenario where the vaccines are allocated during an epidemic and to observe how the allocation stage impacts the optimal allocation strategy and show the dynamic change of the allocation strategy. Thirdly, to evaluate the outcomes of the different strategies, an age-stratified deterministic extended SEIR model was established. The social contact pattern is contained in the model to show how an infectious disease spreads among different age groups. We applied the framework to the COVID-19 vaccine allocation under the different combinations of vaccine parameters and observed some meaningful results; readers can see more details in the section Results and Discussion. Additionally, the framework can be applied to other vaccine allocation situations by changing the evaluation model's parameters or optimizing the objective as their needs.

There are several limitations of this study. First, we investigated vaccine allocation optimization under a fixed social contact pattern. In reality, when an epidemic breaks out, authorities may make policies to let citizens maintain social contact distance, wear masks, and so on to control the spread of the epidemic. Thus, the social contact pattern is dynamic in an epidemic. In addition to age factor, the other factors, such as occupation (42), sex, chronic diseases (16), may relate to susceptibility rate and infected fatality rate, and further influence vaccine allocation. Integrating these factors in this model requires more professional knowledge about COVID-19. These issues should be considered in future work to obtain more realistic results. Second, how to model the effect of a vaccine should be investigated in more detail. In this study, vaccine can decrease the probability of being infected was assumed and vaccines' effect last longer than the time horizon we study, which was simplistic. However, the immunity waning may exist for some vaccines. Third, due to data unavailability, the case study in this study was simulation-based. Otherwise, the test techniques may be immature in some areas at the early stage of vaccine allocation. If these tests are not performed, it is sure that the asymptomatic individuals have opportunity to get vaccine. We leave this scenario as an important consideration for future research.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author/s.

Author contributions

HW and KW contributed to conception and design of the study. HW performed the data analysis and wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

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.934891/full#supplementary-material

Supplementary Figure S1. Pareto optimal fronts when V_e = 20%.

Supplementary Figure S2. Pareto optimal fronts when V_e = 30%.

Supplementary Figure S3. Pareto optimal fronts when V_e = 40%.

Supplementary Figure S4. Pareto optimal fronts when V_e = 50%.

Supplementary Figure S5. Pareto optimal fronts when V_e = 60%.

Supplementary Figure S6. Pareto optimal fronts when V_e = 70%.

Supplementary Figure S7. Pareto optimal fronts when V_e = 80%.

Supplementary Figure S8. Pareto optimal fronts when V_e = 90%.

Supplementary Figure S9. Heatmaps of Pareto optimal solutions when V_e = 20%.

Supplementary Figure S10. Heatmaps of Pareto optimal solutions when V_e = 30%.

Supplementary Figure S11. Heatmaps of Pareto optimal solutions when V_e = 40%.

Supplementary Figure S12. Heatmaps of Pareto optimal solutions when V_e = 50%.

Supplementary Figure S13. Heatmaps of Pareto optimal solutions when V_e = 60%.

Supplementary Figure S14. Heatmaps of Pareto optimal solutions when V_e = 70%.

Supplementary Figure S15. Heatmaps of Pareto optimal solutions when V_e = 80%.

Supplementary Figure S16. Heatmaps of Pareto optimal solutions when V_e = 90%.

Supplementary Figure S17. ICU cases prevention comparation between chosen optimal strategy and baseline.

Supplementary Figure S18. Deaths prevention comparation between chosen optimal strategy and baseline.

Supplementary Figure S19. Pareto optimal fronts of model extension when V_e = 30%.

Supplementary Figure S20. Pareto optimal fronts of model extension when V_e = 60%.

Supplementary Figure S21. Pareto optimal fronts of model extension when V_e = 90%.

Supplementary Figure S22. Heatmaps of Pareto optimal solutions of model extension when V_e = 30%.

Supplementary Figure S23. Heatmaps of Pareto optimal solutions of model extension when V_e = 60%.

Supplementary Figure S24. Heatmaps of Pareto optimal solutions of model extension when V_e = 90%.

References

1. Siettos CI, Russo L. Mathematical modeling of infectious disease dynamics. Virulence. (2013) 4:295–306. doi: 10.4161/viru.24041

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Duijzer LE, van Jaarsveld WL, Wallinga J, Dekker R. Dose-optimal vaccine allocation over multiple populations. Prod Operat Manage. (2018) 27:143–59. doi: 10.1111/poms.12788

PubMed Abstract | CrossRef Full Text | Google Scholar

3. National Academies of Sciences, Engineering, and Medicine. Framework for Equitable Allocation of COVID-19 Vaccine. Washington, DC (2020).

Google Scholar

4. Toner E, Barnill A, Krubiner C, Bernstein J, Privor-Dumm L, Watson M, et al. Interim Framework for COVID-19 Vaccine Allocation and Distribution in the United States. Baltimore, MD: Johns Hopkins Center for Health Security (2020).

5. Andreadakis Z, Kumar A, Román RG, Tollefsen S, Saville M, Mayhew S. The COVID-19 vaccine development landscape. Nat Rev Drug Discov. (2020) 19:305–6. doi: 10.1038/d41573-020-00073-5

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Le TT, Cramer JP, Chen R, Mayhew S. Evolution of the COVID-19 vaccine development landscape. Nat Rev Drug Discov. (2020) 19:667–8. doi: 10.1038/d41573-020-00151-8

PubMed Abstract | CrossRef Full Text | Google Scholar

7. World Health Organization. WHO SAGE Values Framework for the Allocation and Prioritization of COVID-19 Vaccination, 14 September 2020. Geneva: World Health Organization (2020).

Google Scholar

8. Persad G, Peek ME, Emanuel EJ. Fairly prioritizing groups for access to COVID-19 vaccines. JAMA. (2020) 324:1601–2. doi: 10.1001/jama.2020.18513

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Roope LS, Buckell J, Becker F, Candio P, Violato M, Sindelar JL, et al. How should a safe and effective COVID-19 vaccine be allocated? Health economists need to be ready to take the baton. PharmacoEconomics. (2020) 4:557–61. doi: 10.1007/s41669-020-00228-5

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Matrajt L, Longini IM Jr. Optimizing vaccine allocation at different points in time during an epidemic. PLoS ONE. (2010) 5:e13767. doi: 10.1371/journal.pone.0013767

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Medlock J, Galvani AP. Optimizing influenza vaccine distribution. Science. (2009) 325:1705–8. doi: 10.1126/science.1175570

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Preciado VM, Zargham M, Enyioha C, Jadbabaie A, Pappas G. Optimal vaccine allocation to control epidemic outbreaks in arbitrary networks. In: 52nd IEEE Conference on Decision and Control, IEEE (2013).

Google Scholar

13. Yarmand H, Ivy JS, Denton B, Lloyd AL. Optimal two-phase vaccine allocation to geographically different regions under uncertainty. Eur J Operat Res. (2014) 233:208–19. doi: 10.1016/j.ejor.2013.08.027

CrossRef Full Text | Google Scholar

14. Matrajt L, Eaton J, Leung T, Brown ER. Vaccine optimization for COVID-19: who to vaccinate first? Sci Adv. (2020) 7:eabf1374. doi: 10.1101/2020.08.14.20175257

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Tuite A, Fisman DN, Kwong JC, Greer A. Optimal pandemic influenza vaccine allocation strategies for the canadian population. PLoSCurr. (2010) 2:RRN1144. doi: 10.1371/currents.RRN1144

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Mylius SD, Hagenaars TJ, Lugnér AK, Wallinga J. Optimal allocation of pandemic influenza vaccine depends on age, risk and timing. Vaccine. (2008) 26:3742–9. doi: 10.1016/j.vaccine.2008.04.043

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Zhang J, Litvinova M, Liang Y, Wang Y, Wang W, Zhao S, et al. Changes in contact patterns shape the dynamics of the COVID-19 outbreak in China. Science. (2020) 368:1481–6. doi: 10.1126/science.abb8001

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Hao X, Cheng S, Wu D, Wu T, Lin X, Wang C. Reconstruction of the full transmission dynamics of COVID-19 in Wuhan. Nature. (2020) 584:420–4. doi: 10.1038/s41586-020-2554-8

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Ferguson N, Laydon D, Nedjati-Gilani G, Imai N, Ainslie K, Baguelin M, et al. Report 9: impact of non-pharmaceutical interventions (NPIs) to reduce COVID19 mortality and healthcare demand. Imperial College London. (2020) 10:491–7. doi: 10.25561/77482

CrossRef Full Text | Google Scholar

21. Wu JT, Leung K, Leung GM. Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study. Lancet. (2020) 395:689–97. doi: 10.1016/S0140-6736(20)30260-9

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Park M, Cook AR, Lim JT, Sun Y, Dickens BL. A systematic review of COVID-19 epidemiology based on current evidence. J Clin Med. (2020) 9:967. doi: 10.3390/jcm9040967

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Diekmann O, Heesterbeek JAP, Metz JA. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J Math Biol. (1990) 28:365–82. doi: 10.1007/BF00178324

PubMed Abstract | CrossRef Full Text | Google Scholar

24. 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

25. Miettinen K. Introduction to Multiobjective Optimization: Noninteractive Approaches. Multiobjective Optimization. Berlin; Heidelberg: Springer (2008). p. 1–26.

Google Scholar

26. Konak A, Coit DW, Smith AE. Multi-objective optimization using genetic algorithms: a tutorial. Reliability Eng Syst Safety. (2006) 91:992–1007. doi: 10.1016/j.ress.2005.11.018

CrossRef Full Text | Google Scholar

27. Zitzler E, Laumanns M, Thiele L. SPEA2: improving the strength Pareto evolutionary algorithm. TIK-report. (2001). p. 103.

Google Scholar

28. Herrera F, Lozano M, Verdegay JL. Tackling real-coded genetic algorithms: operators and tools for behavioural analysis. Artif Intell Rev. (1998) 12:265–319. doi: 10.1023/A:1006504901164

CrossRef Full Text | Google Scholar

29. Fernandes N. Economic effects of coronavirus outbreak (COVID-19) on the world economy. IESE Business School Working Paper No. WP-1240-E. (2020). doi: 10.2139/ssrn.3557504

CrossRef Full Text | Google Scholar

30. Hou C, Chen J, Zhou Y, Hua L, Yuan J, He S, et al. The effectiveness of quarantine of Wuhan city against the Corona Virus Disease 2019 (COVID-19): a well-mixed SEIR model analysis. J Med Virol. (2020) 92:841–8. doi: 10.1002/jmv.25827

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Pan A, Liu L, Wang C, Guo H, Hao X, Wang Q, et al. Association of public health interventions with the epidemiology of the COVID-19 outbreak in Wuhan, China. JAMA. (2020) 323:1915–23. doi: 10.1001/jama.2020.6130

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Chen S, Chen Q, Yang W, Xue L, Brnighausen T. Buying time for an effective epidemic response: the impact of a public holiday for outbreak control on COVID-19 epidemic spread. Engineering. (2020) 6:1108–14. doi: 10.1016/j.eng.2020.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

33. He S, Peng Y, Sun K. SEIR modeling of the COVID-19 and its dynamics. Nonlinear Dynam. (2020) 101:1667–80. doi: 10.1007/s11071-020-05743-y

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Li R, Chen B, Zhang T, Ren Z, Song Y, Xiao Y, et al. Global COVID-19 pandemic demands joint interventions for the suppression of future waves. Proc Natl Acad Sci. (2020) 117:26151–7. doi: 10.1073/pnas.2012002117

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Chang S, Pierson E, Koh PW, Gerardin J, Redbird B, Grusky D, et al. Mobility network models of COVID-19 explain inequities and inform reopening. Nature. (2021) 589:82–7. doi: 10.1038/s41586-020-2923-3

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Firth JA, Hellewell J, Klepac P, Kissler S, Kucharski AJ, Spurgin LG. Using a real-world network to model localized COVID-19 control strategies. Nat Med. (2020) 26:1616–22. doi: 10.1038/s41591-020-1036-8

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Nicolas H, Martin B, Carlos B, Mark O, Marc M, Sánchez RM, et al. A stochastic agent-based model of the SARS-CoV-2 epidemic in France. Nat Med. (2020) 26:1417–21. doi: 10.1038/s41591-020-1001-6

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Silva PCL, Batista PVC, Lima HS, Alves MA, Guimarães FG, Silva RCP. COVID-ABS: an agent-based model of COVID-19 epidemic to simulate health and economic effects of social distancing interventions. Chaos Solitions Fractals. (2020) 139:110088. doi: 10.1016/j.chaos.2020.110088

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Li R, Pei S, Chen B, Song Y, Zhang T, Yang W, et al. Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2). Science. (2020) 368:489–93. doi: 10.1126/science.abb3221

PubMed Abstract | CrossRef Full Text | Google Scholar

40. China's 2010 Census. Available online at: http://www.stats.gov.cn/tjsj/pcsj/rkpc/6rp/indexch.htm (accessed September 20, 2007).

41. Deb K. Multi-Objective Optimization. Search Methodologies. New York, NY: Springer (2014). p. 403–49. doi: 10.1007/978-1-4614-6940-7_15

CrossRef Full Text

42. Buckner JH, Chowell G, Springborn MR. Dynamic prioritization of COVID-19 vaccines when social distancing is limited for essential workers. Proc Natl Acad Sci. (2021) 118:e2025786118. doi: 10.1073/pnas.2025786118

PubMed Abstract | CrossRef Full Text | Google Scholar

43. 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

44. Chen X, Li M, Simchi-Levi D, Zhao T. Allocation of COVID-19 vaccines under limited supply. MedRxiv. [Preprint]. (2020). doi: 10.1101/2020.08.23.20179820

CrossRef Full Text | Google Scholar

45. Han S, Cai J, Yang J, Zhang J, Wu Q, Zheng W, et al. Time-varying optimization of COVID-19 vaccine prioritization in the context of limited vaccination capacity. Nat Commun. (2021) 12:1–10. doi: 10.1038/s41467-021-24872-5

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Sunohara S, Asakura T, Kimura T, Ozawa S, Oshima S, Yamauchi D, et al. Effective vaccine allocation strategies, balancing economy with infection control against COVID-19 in Japan. PLoS ONE. (2021) 16:e0257107. doi: 10.1371/journal.pone.0257107

PubMed Abstract | CrossRef Full Text | Google Scholar