Optimal vaccine allocation strategy: Theory and application to the early stage of COVID-19 in Japan

: In this paper, we construct an age-structured epidemic model to analyze the optimal vaccine allocation strategy in an epidemic. We focus on two topics: the first one is the optimal vaccination interval between the first and second doses, and the second one is the optimal vaccine allocation ratio between young and elderly people. On the first topic, we show that the optimal interval tends to become longer as the relative efficacy of the first dose to the second dose (RE) increases. On the second topic, we show that the heterogeneity in the age-dependent susceptibility (HS) affects the optimal allocation ratio between young and elderly people, whereas the heterogeneity in the contact frequency among different age groups (HC) tends to affect the effectiveness of the vaccination campaign. A counterfactual simulation suggests that the epidemic wave in the summer of 2021 in Japan could have been greatly mitigated if the optimal vaccine allocation strategy had been taken.


Parameter setting
The unit time is set as 1 day.In Japan, the COVID-19 vaccination program for ordinary people started at April 12, 2021, and the booster (third) vaccination started at December 1, 2021.Because, in this study, we focus on the optimal interval between the first and second doses, we set the time range in our simulation from April 12, 2021 ( 0) to November 30, 2021 (  232).As the vaccine age does not exceed the calendar time, we set  232.Let the unit age be 1 year, and the age interval be 0,100 , that is,  100.We set  0.2 and  0.1 so that the average incubation period is 1/ 5 days [1,2] and the average infectious period is 1/ 10 days [3].We assume that the severe class is composed of individuals who will die due to COVID-19, and set   as shown in Figure S2 by dividing the cumulative deaths by the cumulative cases in each age class as of November 30, 2021, using the open data in [4].The vaccine hesitancy rate ℎ  is estimated as in Figure S3 using the vaccination data in [5].
For the initial condition, let be the population age distribution in Japan as of April 2021 ( 0).Using the data in [6], we fix   as shown in Figure S4.Here, we normalize   so that    1.Each population then describes a proportion to the total population.For simplicity, we assume that   0, and fix   ,   and   as shown in Figure S5 using the data in [4].  can then be calculated as         .We assume that the infection rate  ,  has the following form: We can numerically compute   by using a discretization method as in [8].If we consider the case of ℛ 1.5, then  is modified so that   1.5 is attained (note that   is proportional to ).
Let  be the length of the vaccination interval between the first and second doses.Assuming that severe individuals would not be newly vaccinated, the number of the first vaccination shots at time  is calculated as  ,   ,   ,   ,   ,   , where  1.26 10 is the total population in Japan as of 2021 [6].The number of the second vaccination shots at time  is given by  , ,   , ,   , ,   , ,   .
We assume that the vaccination rate is separable: (1)  ,      .
It then follows that where  and  ( 1,2) denote the efficacy of the -th vaccination in reducing the risk of infection and the risk of disease-related death, respectively.Here, according to [9], we assume that the vaccine efficacy linearly decreases with waning rate  0. In this study, we fix  1/600 so that the efficacy decreases to its half after 300 days passed [9].We consider the Pfizer and AstraZeneca vaccines taking the mean of a dataset in [10]: ,  ,  ,  0.63, 0.90, 0.80, 0.94 , Pfizer, 0.62, 0.64, 0.80, 0.85 , AstraZeneca.

Objective functions
To evaluate the effectiveness of vaccination program, we calculate the total number of cases as To compute the reduction ratio, we divide each of these functions by those without vaccination (  0).

Counterfactual simulation for Japan in 2021
For the baseline scenario, we assume that the vaccination interval between the first and second doses is 3 weeks, the vaccine efficacy is as of Pfizer, the ratio  of vaccine allocation to the elderly people is 0.9, and both  and  are medium.We then estimate the time-varying infection rate with    by fitting the following function to the daily reported cases in Japan: where  is the detection ratio, that is, the ratio at which a newly infected individual is finally reported.More precisely, to estimate the infection rate on a day, we compare the simulation result and the real data of the newly reported cases in past 7 days, and find the parameter that minimizes the sum of the squares error.In our simulation, we assume that  0.5.The comparison of the curve of   and the real data [4] is shown in Figure S9, top.Note that this function   is uniquely calculated.For this estimated infection rate, we compare the daily number of newly reported deaths     ,   ,   ,   , and the real data [4] as shown in Figure S9, bottom.Here, for a technical reason on the fitting of   , we compare the total deaths in the period from May 15, 2021 to November 30, 2021.
Regarding this setting as a baseline, we investigate how the total deaths could be reduced if the optimal vaccination interval and the optimal ratio of allocation to the elderly people had been taken.

Figure S1 .
Figure S1.Transfer diagram of our model.

Figure S2 .
Figure S2.Age-specific probability at which a recovered individual becomes severe.

Figure S3 .
Figure S3.Age-specific vaccine hesitancy rate h a .

Figure S4 .
Figure S4.Initial population age distribution P a .

Figure S5 .
Figure S5.Initial age distributions of exposed   , infectious   and recovered   populations.

Figure S6 .
Figure S6.Age dependent susceptibility β a (HS) and contact frequency β x among individuals whose age difference is x (HC).Three levels of heterogeneity (high, medium and low) are considered to each of them.

Figure S7 .
Figure S7.Daily number of vaccination shots in Japan from April 12, 2021 to November 30, 2021.

Figure S9 .
Figure S9.Daily number of newly reported cases (top) and deaths (bottom).Functions   and   are derived from our model.

Figure S10 .
Figure S10.The main part of our simulation code for MATLAB.Some symbols are defined in other parts or files.

Table S1 .
The parameter values in our model.
[7]re  ,  denotes the probability density function of normal distribution with mean  and standard derivation .represents the heterogeneity in the age-dependent susceptibility (HS) and   represents the heterogeneity in the contact frequency among different age classes.The reason why we assume that   is higher in younger age group is that elderly people might be more careful and more likely to reduce the contact opportunity because the mortality of COVID-19 is quite high in those people.The reason why we assume that   has five peaks is that the contact opportunity among children and their parents and grandparents might be high.In addition, we assume that if   (HC) is high, then people become more likely to contact with people in similar age group, whereas if   (HC) is low, then people become more likely to contact beyond age groups and the mixing becomes more homogeneous.The parameter  is modified to fix the basic reproduction number ℛ for different   and   .Following the classical theory[7], ℛ is defined by the spectral radius   of the following next generation operator : ,       ,where  is the infection transmission rate,   is the age-dependent susceptibility, and   is a distance function representing the contact frequency among individuals whose age difference is .More precisely, according to the three levels of heterogeneity (high, medium and low), we set   and   as follows (see FigureS6): [5]re   denotes the total number of vaccination at time .We fix   as shown in FigureS7using the data of vaccine distribution for COVID-19 in Japan from April 12, 2021 to November 30, 2021[5].  can be used to incorporate the priority of vaccine allocation to individuals aged .
Let ℓ  be the average life expectancy at age , which is estimated from the data in[11]as shown in FigureS8.The total number of disease-induced deaths weighted by the average life expectancy is as follows: