Evaluating the COVID-19 vaccination program in Japan, 2021 using … – Nature.com

Conversion to infections

COVID-19 was designated a notifiable disease under the infectious disease law of Japan as of 2021. All individuals suspected of being infected with SARS-CoV-2 were tested via PCR or quantitative antigen test at medical facilities. They were then requested to remain in home isolation and undergo investigation by municipal public health centers to identify their close contacts. Information of confirmed cases (e.g., age and sex) was registered in the Health Center Real-time Information-sharing System on COVID-19 (HER-SYS) by medical facilities or municipal public health centers. Supplementary Fig. S1A shows the number of confirmed cases from the beginning of the primary series (the first and second doses) of the vaccination program through the end of November 2021. In the end of November 2021, SARS-CoV-2 in Japan was dominated by Delta variant to which the vaccine effectiveness was known to have been greatly diminished, sometimes by 10%, compared with other variants that circulated earlier47,48,49.

The time of infection for all confirmed COVID-19 cases retrieved from HER-SYS was backcalculated using a previously estimated distribution of the interval between infection and illness onset, assumed to follow a log-normal distribution with a mean of 4.6days and standard deviation (SD) of 1.8days50, 51. Cases without a date of symptom onset were backcalculated using the time difference from symptom onset to reporting, assumed to follow a log-normal distribution with a mean of 2.6days and SD of 2.1days, as previously estimated using cases with information for the date of symptom onset. Non-parametric backcalculation was performed using the R-package surveillance (version 1.20.3). To address the issue of reporting bias, we explored different reporting coverages: 0.125, 0.25, 0.5, and 1.0 (no bias) by multiplying the backcalculated cases by 1 and dividing by reporting coverage to finally obtain the number of infections.

SARS-CoV-2, all vaccinated individuals retrieved from the Vaccine Record System (VRS) were converted into immunized people according to time. The data comprised the sex, age, and date of vaccination for vaccinated individuals. We assumed that all people who received the first dose were subsequently vaccinated with the second dose at an interval of 21days (Supplementary Fig. S1B). According to statistics of the VRS, there was a very small discrepancy in vaccination coverage between the first dose (75.19%) and the second dose (74.61%) as of the end of December 202152; therefore, we could obtain a certain consensus on the usage data for people vaccinated with the first dose only. For the conversion, we used a profile of vaccine efficacy involving waning immunity for the primary series used by Gavish et al.19, which was based on previous estimates53,54. Given the widespread use of the messenger RNA vaccine BNT162b2 (Pfizer/BioNTech) in Japan (more than 80% of individuals received this vaccine by the end of November 2021)23, we assumed that published estimates could directly be applied to the case of Japan. Further details and background of the primary series in Japans vaccination program are described elsewhere17.

To adapt the following transmission model, we used the number of vaccinated individuals and the profile of vaccine efficacy to estimate the immune fraction in age group (a) at calendar time (t), ({l}_{a,t}), which is expressed as:

$${l}_{a,t}=frac{1}{{n}_{a}}sum_{s=1}^{t-1}{v}_{a,t-s}{h}_{s}$$

(1)

where ({n}_{a}) is the population size in age group (a) in 202155, ({v}_{a,t}) denotes the number of vaccinated individuals in age group (a) at calendar time (t), and ({h}_{s}) represents the vaccine profile. Supplementary Fig. S2 displays the estimated immune fraction by age group.

We developed the time-dependent transmission model that accounts for heterogeneous transmission between age groups, fitting the model to observed incidence data and estimating unknown parameters. We used the following renewal equation to infer the transmission dynamics underlying the COVID-19 epidemic, which is described as:

$${i}_{a,t}=sum_{b=1}^{10}sum_{tau =1}^{t-1}{{varvec{R}}}_{{varvec{a}}{varvec{b}},{varvec{t}}}{i}_{b,t-tau }{g}_{tau },$$

(2)

where ({i}_{a,t}) represents the number of infections with SARS-CoV-2 in age group (a) at day (t) and ({g}_{tau }) indicates the probability density function of the generation interval, assumed to follow a Weibull distribution with a mean of 4.8days and SD of 2.2days51,56. ({{varvec{R}}}_{{varvec{a}}{varvec{b}},{varvec{t}}}) denotes the effective reproduction number, interpreted as the average number of secondary cases in age group (a) generated by a single primary case in age group (b) at calendar time (t). To capture the impact of vaccination, ({{varvec{R}}}_{{varvec{a}}{varvec{b}},{varvec{t}}}) was decomposed as:

$${{varvec{R}}}_{{varvec{a}}{varvec{b}},{varvec{t}}}=left(1-{l}_{a,t}-frac{sum_{k=1}^{t-1}{i}_{a,k}}{{n}_{a}}right){{varvec{K}}}_{{varvec{a}}{varvec{b}}}p{h}_{t}{d}_{t}{c}_{t}$$

(3)

where (sum_{k=1}^{t-1}{i}_{a,k}) represents the cumulative number of previous infections after 16 February 2021. ({{varvec{K}}}_{{varvec{a}}{varvec{b}}}) is considered a next-generation matrix, which was modeled as ({{varvec{K}}}_{{varvec{a}}{varvec{b}}}={{{s}}}_{{{a}}}{{varvec{m}}}_{{varvec{a}}{varvec{b}}}), where ({{{s}}}_{{{a}}}) represents relative susceptibility and ({{varvec{m}}}_{{varvec{a}}{varvec{b}}}) denotes the contact matrix; we rescaled a previously quantified next-generation matrix during the initial phase of the COVID-19 epidemic in 2021 attributable to the Alpha variant57. Because the oldest age group was65years in the previous estimate, we reconstructed the epidemic curve with new age groups: 09, 1019, 2029, 3039, 4049, 5059, 6069, 7079, 8089 and90years and estimated ({{varvec{K}}}_{{varvec{a}}{varvec{b}}}) by fitting the model to observed cases (Supplementary Fig. S3). The detailed methods are explained elsewhere57,58. We assumed that the contact rates among groups aged70years were the same as those aged65 in the contact matrix, ({{varvec{m}}}_{{varvec{a}}{varvec{b}}}), which is based on a social epidemiological survey conducted prior to the COVID-19 pandemic in Japan42. With respect to the above explanation, those early terms in Eq.(3) could capture the effective heterogeneous interactions between infectees and infectors, which accounts for the immune fraction owing to vaccination and infections among susceptible individuals (i.e., infectees). (p) denotes the scaling parameter involving all terms in Eq.(3) and ({h}_{t}) expresses the change in mobility. The variable, ({h}_{t}), related to human mobility was decomposed as:

$${h}_{t}={omega }^{community}{alpha }_{t}^{community}+{omega }^{house}{alpha }_{t}^{house}+{omega }^{work}{alpha }_{t}^{work},$$

(4)

where (omega) means the coefficient of human mobility in the community, household, or workplace relative to the community setting (i.e., ({omega }^{community}) is equivalent to 1). The coefficient, ({alpha }_{t}), describes a proxy of the intensified contacts in three different settings retrieved from Google's COVID-19 community mobility report in Japan59. Those data were smoothed using a 7-day moving average (Supplementary Fig. S4). ({d}_{t}) represents the increase in transmissibility of the Delta variant compared with earlier variants, which was formulated as ({d}_{t}=r{u}_{t}), where (r) is the scaling parameter for transmissibility and ({u}_{t}) represents the profile of increased transmissibility. We assumed that ({u}_{t}) increased with the detected proportion of COVID-19 cases owing to the Delta variant in Japan60, which was modeled using a logistic curve. We then rescaled ({u}_{t}) up from 1 to a maximum of 1.561,62,63. A comparison between the predicted and observed proportion is shown in Supplementary Fig. S5. ({d}_{t}) was parameterized as 1 before 20 May 2021, when we assumed that the proportion of infections with the Delta variant started to increase at population level. Finally, ({c}_{t}) expresses the influence of consecutive holidays, defined as more than 3days in the present study. Moreover, we added Obon season, the national religious season associated with Buddhist tradition, to those holidays. Not all consecutive days in this period (from 13 to 16 August 2021) were regarded as holidays; however, many Japanese people travel and/or visit their relatives during this season. We modeled ({c}_{t}) as ({c}_{t}=e{beta }_{t}), where (e) accounts for the coefficient of holiday influence and ({beta }_{t}) was assigned 1 if the day was aligned with consecutive holidays; otherwise ({c}_{t}) was parameterized as 1.

We first computed the counterfactual scenario, i.e., without vaccination. We also explored three additional hypothetical scenarios: (1) the vaccination program was implemented sooner than the actual program, reaching a maximum number of vaccinated individuals 14days earlier than the observed pace (hereafter early schedule scenario); (2) the vaccination schedule was delayed, reaching a peak in the number of vaccinated people 14days slower than the observed pace (late schedule scenario); and (3) adolescents and people aged 1059years were vaccinated more and faster (elevated scenario). To explore different counterfactual scenarios, we first regressed the vaccination coverage using the logistic function by age group, which is modeled as:

$$E({v}_{a,t})=frac{{pi }^{1}}{1+mathrm{exp}(-{pi }^{2}(t-{pi }^{3}))},$$

(7)

where ({pi }^{1}), ({pi }^{2}), and ({pi }^{3}) represent the carrying capacity (eventual coverage of the primary series), speed of increase in the vaccination coverage, and requisite duration for the half coverage of ({pi }^{1}) (also representing the peak day for the number of vaccinated individuals), respectively. We performed maximum likelihood estimation to estimate ({pi }^{1}), ({pi }^{2}), and ({pi }^{3}) by age group. Comparisons between the predicted and observed number of vaccinated people by age group are shown in Supplementary Fig. S6.

We assumed that the days with the maximum number of vaccinated people (i.e., days that 50% of the carrying capacity was achieved) were 14days earlier in the Early scenario and later in the Late scenario than the observed. For the Elevated scenario, we assumed that people aged 1059years had earlier peaks in the number of vaccinated individuals, as with the Early scenario. Additionally, people aged 1019years and aged 2049years were assumed to reach 70% and 90% in eventual vaccination coverage (({pi }^{1})), respectively. People aged50years had already reached more than 90% of the vaccination coverage by the end of November 2021. We did not consider vaccination among individuals aged less than 10years because children were not eligible to be vaccinated during the primary series of the program in Japan. All scenarios of the vaccination program by age group are shown in Supplementary Fig. S7.

We assumed that the daily counts of infections followed a Poisson distribution, and the likelihood function with unknown parameters, (theta ={p,{omega }^{house},{omega }^{work},r,e}), was represented as:

$$Lleft(theta ;{i}_{a,t}right)=prod_{t}prod_{a}frac{{E({i}_{a,t})}^{{i}_{a,t}}mathrm{exp}(-E({i}_{a,t}))}{{i}_{a,t}!}$$

(8)

By minimizing the loglikelihood function, we estimated (theta). The 95% confidence intervals (CI) were calculated from 1000 bootstrap iterations using the multivariate normal distributions of the parameters. We estimated a series of parameters by reporting coverage in the present study. All estimated parameters with 95% CIs are shown in Supplementary Table S1. Supplementary Fig. S8 demonstrates the fitting outcome of the predicted and observed infections with SARS-CoV-2 by age group, with reporting coverage of 1 (i.e., no ascertainment bias). Supplementary Fig. S9 compares the predicted and observed infections by reporting coverage.

Using the estimated parameters, (theta), we explored hypothetical scenarios by varying the timing and the recipients of vaccination. For this, we used infections already backcalculated 14days back from the start of vaccination as the initial condition.

Because the effective reproduction number in Japan conventionally uses an estimate for the entire population, we also calculated an effective reproduction number based on the total number of cases at calendar time (t), ({R}_{t}), in each counterfactual scenario using the total number of infections with SARS-CoV-2. Using an equation similar to Eq.(2), the total number of infections, ({i}_{t}^{total}), was modeled as:

$${i}_{t}^{total}={R}_{t}sum_{tau }^{1-tau }{i}_{t-tau }^{total}{g}_{tau }.$$

(9)

Assuming the daily case counts followed a Poisson distribution, we estimated ({R}_{t}) using maximum likelihood estimation51.

To compute the mortality impact, we estimated the age-specific infection fatality risk (IFR) according to reporting coverage in the present study. First, we formulated the cumulative number of deaths in age group (a) resulting from cases infected during the research period in unvaccinated and vaccinated individuals, respectively, which are described as:

$$left{begin{array}{c}{D}_{a}^{unvaccinated}={IFR}_{a}sum_{t=17mathrm{ Feb }2021}^{30mathrm{ Nov }2021}(1-{epsilon }_{a,t}){widehat{i}}_{a,t}\ {D}_{a}^{vaccinated}={IFR}_{a}(1-{VE}_{a})sum_{t=17mathrm{ Feb }2021}^{30mathrm{ Nov }2021}{epsilon }_{a,t}{widehat{i}}_{a,t}end{array}right.,$$

(10)

where ({epsilon }_{a,t}) represents the time-varying proportion of vaccinated people among confirmed cases in age group (a) at calendar time (t), and ({widehat{i}}_{a,t}) is the expected number of infections estimated from the transmission model. ({VE}_{a}) expresses the vaccine-induced reduction in mortality estimated in 2021 for Japan64. We obtained ({epsilon }_{a,t}) by modeling cases with a vaccination history registered in HER-SYS using a logistic function. The observed proportion was calculated as 7-day moving average and shifted5days because of the conversion for the time of infection. Also, to account for the age groups used in the present study, people aged 1019, 2029, 3039, 4049, 5059 and60years were utilized as people aged 1524, 2534, 3544, 4554, 5564 and65years for the proportion retrieved from HER-SYS, respectively. Supplementary Fig. S10 shows the comparison between the model prediction and observed proportions.

To estimate IFR by age group, the following likelihood equation was used:

$$Lleft({lambda }_{a};{widehat{i}}_{a,t},{D}_{a}right)=prod left(begin{array}{c}sum {widehat{i}}_{a,t}\ {D}_{a}end{array}right){lambda }_{a}^{{D}_{a}}{left(1-{lambda }_{a}^{{D}_{a}}right)}^{sum {widehat{i}}_{a,t}{-D}_{a}},$$

(11)

where ({lambda }_{a}) denotes the risk of death in age group (a), modeled as:

$${lambda }_{a}=frac{({D}_{a}^{unvaccinated}+{D}_{a}^{vaccinated})}{sum {widehat{i}}_{a,t}}.$$

(12)

({D}_{a}) is the cumulative number of deaths reported from 10 March to 21 December 2021 in age group (a), which was retrieved from the Ministry of Health, Labour and Welfare of Japan, accounting for the reporting delay of 21days.65 By minimizing the negative logarithm of Eq.(8), we estimated ({IFR}_{a}). We performed this process for each reporting coverage. Supplementary Fig. S11 displays the estimated IFR by reporting coverage and age group. Finally, we estimated the cumulative number of deaths as an aggregation of ({D}_{a}^{unvaccinated}) and ({D}_{a}^{vaccinated}) in Eq.(7) according to different counterfactual scenarios of varying ({widehat{i}}_{a,t}). We only applied the first equation in Eq.(7), i.e., ({D}_{a}^{unvaccinated}), for the counterfactual scenario in the absence of vaccination.

For calculation of the death toll, we altered only a parameter representing the requisite duration for the half coverage of a carrying capacity to coincide with changes in vaccine recipients in the counterfactual vaccination scenarios. Because of this exercise, we were able to model the specific proportion of vaccinated people among confirmed cases according to different vaccination scenarios. The principal idea of the logistic model is explained in the early subsection.

This study was conducted according to the principles of the Declaration of Helsinki. Informed consent was obtained for reporting the diagnosis. The authors did not have an access to any individual identity information, and this research was approved by the Ethics Committee of Kyoto University Graduate School of Medicine (approval number R2673).

Read more from the original source:

Evaluating the COVID-19 vaccination program in Japan, 2021 using ... - Nature.com