Modeling and predicting individual variation in COVID-19 vaccine-elicited antibody response in the general population

Abstract

As we learned during the COVID-19 pandemic, vaccines are one of the most important tools in infectious disease control. To date, an unprecedentedly large volume of high-quality data on COVID-19 vaccinations have been accumulated. For preparedness in future pandemics beyond COVID-19, these valuable datasets should be analyzed to best shape an effective vaccination strategy. We are collecting longitudinal data from a community-based cohort in Fukushima, Japan, that consists of 2,407 individuals who underwent serum sampling two or three times after a two-dose vaccination with either BNT162b2 or mRNA-1273. Using the individually reconstructed time courses of the vaccine-elicited antibody response based on mathematical modeling, we first identified basic demographic and health information that contributed to the main features of the antibody dynamics, i.e., the peak, the duration, and the area under the curve. We showed that these three features of antibody dynamics were partially explained by underlying medical conditions, adverse reactions to vaccinations, and medications, consistent with the findings of previous studies. We then applied to these factors a recently proposed computational method to optimally fit an “antibody score”, which resulted in an integer-based score that can be used as a basis for identifying individuals with higher or lower antibody titers from basic demographic and health information. The score can be easily calculated by individuals themselves or by medical practitioners. Although the sensitivity of this score is currently not very high, in the future, as more data become available, it has the potential to identify vulnerable populations and encourage them to get booster vaccinations. Our mathematical model can be extended to any kind of vaccination and therefore can form a basis for policy decisions regarding the distribution of booster vaccines to strengthen immunity in future pandemics.

Data Availability:The research ethics review application form submitted to and approved by the Ethics Committee of Hirata Central Hospital (number 2021-0611-1) and Fukushima Medical University School of Medicine (number 2021-116) states that the anonymized clinical dataset can only be accessed by the research groups of this study, as it contains potentially identifying and sensitive patient information. In accordance with the PLOS Data Policy, which provides the option of de-identifying human research participant data and adding noise to protect participant privacy, we provide a synthetic patient dataset created using CTGAN (https://github.com/sdv-dev/CTGAN) based on the original dataset (S1 Data). This synthetic patient dataset, while not containing personally identifiable information, retains the same statistical properties as the original and can be analyzed in the same manner as the original.

Funding:This study was supported by Medical & Biological Laboratories Co., Ltd., Shenzhen YHLO Biotech Co., Ltd., the distributor and manufacturer of the antibody measurement system (iFlash 3000), Research Center for Advanced Science and Technology in the University of Tokyo, and in part by Scientific Research (KAKENHI) B 18H01139 (to S.I.), 16H04845 (to S.I.); Scientific Research in Innovative Areas 20H05042 (to S.I.); AMED Strategic International Brain Science Research Promotion Program 22wm0425011s0302 (to K.A.); AMED JP22dm0307009 (to K.A.); AMED CREST 19gm1310002 (to S.I.); AMED Development of Vaccines for the Novel Coronavirus Disease, 21nf0101638h0001 (to M.T.), 21nf0101638s0201 (to S.I.); AMED Japan Program for Infectious Diseases Research and Infrastructure, 22wm0325007h0001 (to S.I.), 22wm0325004s0201 (to S.I.), 22wm0325012s0301 (to S.I.), 22wm0325015s0301 (to S.I.); AMED Research Program on Emerging and Re-emerging Infectious Diseases 22fk0108140s0802 (to S.I.); AMED Research Program on HIV/AIDS 22fk0410052s0401 (to S.I.); AMED Program for Basic and Clinical Research on Hepatitis 22fk0210094 (to S.I.); AMED Program on the Innovative Development and the Application of New Drugs for Hepatitis B 22fk0310504h0501 (to S.I.); AMED Research Program on Hepatitis 19fk0210036h0502 (to S.I.), 19fk0310114h0103 (to S.I.); JST MIRAI JPMJMI22G1 (to S.I.); Moonshot R&D JPMJMS2021 (to K.A. and S.I.) and JPMJMS2025 (to S.I.); The National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (2022R1C1C2003637) (to K.S.K.); The Shinnihon Foundation of Advanced Medical Treatment Research (to S.I.); SECOM Science and Technology Foundation (to S.I.); The Japan Prize Foundation (to S.I.); and Kowa Company, Ltd. (to M.T.). URL of funder websites: Medical & Biological Laboratories Co., Ltd. https://www.mblbio.com/e/ Shenzhen YHLO Biotech Co., Ltd. https://en.szyhlo.com/ Research Center for Advanced Science and Technology in the University of Tokyo https://www.rcast.u-tokyo.ac.jp/en/ KAKENHI https://www.jsps.go.jp/english/e-grants/ AMED https://www.amed.go.jp/en/index.html Moonshot R&D https://www.jst.go.jp/moonshot/en/ The National Research Foundation of Korea https://www.nrf.re.kr/eng/main The Shinnihon Foundation of Advanced Medical Treatment Research http://www.shinnihon-zaidan.jp/ The Japan Prize Foundation https://www.japanprize.jp/en/index.html Kowa Company, Ltd. https://www.kowa.co.jp/

Competing interests: YKaneko is employed by Medical & Biological Laboratories, Co. (MBL, Tokyo, Japan). MBL imported the testing material used in this research. YKaneko participated in the testing process; however, he did not engage in the research design and analysis. YKobashi and MT received a research grant from Pfizer Health Research Foundation for research not associated with this work.

Introduction

The global distribution of vaccines for coronavirus disease 2019 (COVID-19) and the high vaccine potency and coverage will bring the pandemic caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) under control. As we learned from the COVID-19 pandemic, vaccination is an important part of a multi-faceted public health response to future pandemic illnesses. An unprecedentedly large volume of high-quality data have been accumulated during the COVID-19 pandemic, and we should use these data to prepare for future pandemics. In particular, given the limited global supply of vaccines during the early phase of a pandemic, determining vaccination priority is important to effectively distribute doses and achieve early disease control. In addition, waning of vaccine efficacy becomes a major concern during vaccination campaigns, as we observed for COVID-19 vaccinations [1,2]. Because a rapid decline in vaccine-elicited antibodies may result in breakthrough infections [3,4], additional, appropriately timed booster vaccinations will be required to maintain a high-level antibody response at both individual and population levels. Thus, to shape an effective vaccination strategy, it is important to quantify the individual-level time course of antibody dynamics and to identify “vulnerable” populations with sustained low antibody titers. Thus, by using a large volume of high-quality data on COVID-19 vaccination, our purpose in this study was to establish an approach to quantifying vaccine-elicited time-course antibody dynamics and predicting individual-level antibody responses with a simple and noninvasive method.

To date, the findings of studies evaluating vulnerable populations for COVID-19 have been biased by the selection of study population. Rather than studying a general population, studies to date have focused on specific subpopulations stratified by variables like age, sex, lifestyle habits, comorbidities, adverse reactions, or medication use [2,3,5–8]. However, variation in vaccine-elicited antibody responses among healthy individuals per se has also been suggested to correlate with breakthrough infection risk [2,9–11]. In fact, although the heterogeneity in antibody responses over time is important for identifying the characteristics of vulnerable populations (in addition to standard risk factors such as age and comorbidities) [2,3,5–8], individual-level variation remains poorly understood. This is also true for viruses beyond SARS-CoV-2, because the dynamics of antibody responses in humans has not been described in detail for any vaccination.

Here we used a mathematical model to describe the process of differentiation from naïve B cells to plasma cells to accurately reconstruct individual vaccine-elicited antibody dynamics. To fit the model, we used longitudinal antibody measurements from non-sequential and sequential blood sampling in the Fukushima vaccination cohort (a community-based cohort in Fukushima, Japan). Of note, because our vaccination cohort consisted of participants from a primarily rural area where the prevalence of COVID-19 was relatively low, the samples from this cohort are ideal for modeling vaccine-elicited antibody dynamics. Thanks to the reconstructed time-course of antibody dynamics, we were able to compare the antibody response at the same time points at an individual level. This overcomes another limitation of current vaccination studies, which can only directly compare antibody titers at a population level rather than at an individual level at specific time points (i.e., date on blood sampling) [2,3,5]. The model parameters describe highly variable individual-level antibody responses, allowing us to partially predict variation in vaccine response on the basis of personal information including age, adverse reactions, comorbidities, and medication use. Furthermore, we devised a personalized antibody score that aims to identify individuals with higher or lower antibody titers from their personal information. Although this score may not be able to detect all individuals, it has the potential to be used by medical practitioners to encourage individuals with low predicted antibody levels to get booster vaccinations. We stress that our approach will be easily applied to reconstruct antibody responses even after the third, fourth and fifth booster doses.

Results

The Fukushima vaccination cohort

Our vaccination cohort, the Fukushima vaccination cohort, was conducted beginning in April 2021 and consisted of participants from a primarily rural area where COVID-19 prevalence was relatively low: Soma City, Minami Soma City, and Hirata village in Fukushima in Japan (Fig 1A). The data used in this study were obtained from April 2021 through December 2021. The participants included health care workers, frontline workers, administrative officers, general residents, and residents of long-term care facilities. In total, 2,526 participants who had been vaccinated with the Pfizer BNT162b2 or Moderna mRNA-1273 vaccine were recruited, and 2,407 participants were included in the final data analysis (see Fig 1B and Methods for more details). The age and sex distributions of the participants are shown in Fig 1C, and the sample characteristics and information on adverse vaccine events stratified by age (with p-values) are provided in Table 1. A portion of this cohort was described previously for the time period extending to 6 months after the first dose of mRNA vaccine [12–19].

Fig 1. Characteristics of the Fukushima vaccination cohort: (A) Locations of Soma city, Minami Soma city, and Hirata village in Fukushima prefecture are described. (B) Flowchart of the vaccination cohort along with the number of participants, testing, and inclusion criteria for our analysis are described. (C) Age and sex distributions in the cohort are shown. (D) Timeline of sample collection for each cohort participant (N = 2,526 participants, 5195 total samples) is described. The timings of blood samplings are indicated in black circles. Shaded areas indicate early ([less than]89 days), middle (90–179 days), and late (>180 days) time periods after the first vaccination. Dates for the second vaccination are shown as the distributions in the bottom panel. (E) Vaccination and blood sampling periods in the cohort along with the number of COVID-19 cases (i.e., cases) in Soma city, Minami Soma city, and Hirata village [67,68] are shown. (F) Longitudinal IgG(S) and neutralizing activity measured by CLIA are separately plotted by time after the first vaccination and age.

Table 1. Basic demographics for the Fukushima vaccination cohort.

Here we investigated antibody titers in the Fukushima vaccination cohort in individuals sampled longitudinally (serum was collected at 2 or 3 different timepoints) for around 4 to 9 months after the second primary dose of mRNA vaccine (see Fig 1DE for details). Notably, the number of SARS-CoV-2 infections in this rural area was extremely low (Fig 1E), so that we could minimize the influence of natural breakthrough infections. Compared with some of the largest cohorts in the world [1,5,20–23], the Fukushima vaccination cohort is community-based, includes non-health workers, has very few dropouts among more than 2,000 individuals who were consecutively sampled (only 3.3%), includes all necessary information for all participants, and includes measures of several modalities of antibody titers including neutralizing activity.

We performed chemiluminescent immunoassay (CLIA) to measure antibody titers as a measure of humoral immune status after the first COVID-19 vaccination (i.e., a total of 5195 IgG(S), 5195 neutralization activity, and 4969 IgG(N) assays were performed) (Methods). Fig 1F shows the overall profile of IgG(S) and neutralization activity against the Wuhan strain in this study. We investigated longitudinal data for IgG(S) in the same individuals because IgG(S) covers a wider range of antibody responses and is more sensitive than neutralization activity. In fact, these two measurements are highly correlated with each other (correlation coefficient of 0.93) (S1 Fig), and previous studies showed that neutralizing antibody and IgG(S) titers correlate with vaccine-mediated protection, even against variants of concern (i.e., vaccine efficacy) [9,24,25]. This is because vaccines containing the original Wuhan virus spike protein induce variant-reactive memory B cells targeting multiple variants of concern, including the Omicron variant [26]. To prepare for a rapid response to an early phase of a future pandemic, we hereafter investigated longitudinal data for IgG(S) titers against the ancestral strain in the same individuals, which we used as a biomarker for vaccine-elicited immune response at the beginning of a COVID-19 vaccination program, and we used these data to develop an approach for establishing an "antibody score".

Deriving measures of peak, duration, and area under the curve of vaccine-elicited antibody dynamics

We developed a mathematical model describing the vaccine-elicited antibody dynamics to evaluate the impact of primary two-dose COVID-19 vaccination on rapid immunity at the individual level. We fully reconstructed the dynamics of IgG(S) titers after the first vaccination for 2,407 individuals in the Fukushima vaccination cohort in S2 Fig (see Methods in detail). Of note, we validated our mathematical model and the reconstructed antibody titer curves using independent datasets (see Methods for details). Then we extracted the “features” described in Fig 2A for each individual: the peak, duration, and area under the curve (AUC) of the reconstructed antibody dynamics. To quantify these features, we here assumed ATH = 100, and determined ts and te corresponding to the time for the antibody titer to be greater than and smaller than ATH, respectively. ts and te are calculated as the minimum and maximum values of the time, respectively, during which the reconstructed IgG(S) remains above ATH. Therefore, the duration and AUC of the antibody titer are formulated by te − ts and , respectively. In addition, defining tp to be the time for the antibody titer to reach its peak, the peak titer is A(tp). For 89 individuals whose peak titer was below ATH, the duration and AUC were both determined to be 0. In Fig 2B, we summarized distributions of the AUC, duration, and peak for 2,407 participants. Note that a similar trend was obtained under different ATH. These features allowed us to quantitatively compare vaccine-elicited antibody dynamics among the participants (see next section).

Fig 2. Quantifying vaccine-elicited antibody dynamics: (A) Vaccine-elicited antibody response after the first vaccination (i.e., t = 0) is described with the following “features”: the peak (A(tp)), duration (te − ts), and AUC of the antibody titers. The vertical and horizontal dashed lines correspond to the date of the second vaccination and the arbitrary threshold (ATH) for calculating the duration and AUC, respectively. (B) Distributions of the extracted features from the reconstructed antibody dynamics (i.e., the peak, duration, and AUC) for 2,407 participants are plotted. The dataset for each distribution was normalized by the value corresponding to the 95th percentile of data values, and data with values larger than this value were removed to improve the visibility of the figure.

Characterizing vaccine-elicited antibody dynamics

To see how individual background factors contributed to the three features, we divided the peak-duration plane into four quadrants (groups 1–4) by taking the median values of peak and duration of titer as cutoffs (Fig 3A). We collected basic demographic and health information from the participants, including underlying medical conditions, adverse reactions to vaccinations, and medications, as described in Table 1. We then investigated whether there were differences in these basic variables among the four groups.

Fig 3. Characterizing and scoring vaccine-elicited antibody dynamics: (A) The peak-duration plane was divided into four quadrants (groups 1 to 4). The cutoff values for the peak (log10) and the duration were 2.94 and 208, respectively. The distribution of 2,407 participants on the plane is shown as a density plot. (B) Distributions of age, BMI and the interval between the two vaccine doses in groups 1 to 4 are plotted. (C) The frequency of each categorical variable in groups 1 to 4, stratified by gender, is shown as a heatmap. In the figure, 1, 2, 3, 4 refer to groups 1–4 in men, whereas 5, 6, 7, 8 refer to groups 1–4 in women.

We first analyzed continuous variables such as age, BMI, and the interval between the two vaccinations. We found that group 1 (high peak, long duration) had the youngest participants on average, while group 4 (low peak, short duration) had the oldest participants on average (Fig 3B top). Groups 1 and 3 (high peak) had longer intervals between the two vaccinations than groups 2 and 4 (low peak) (Fig 3B middle). In contrast, there was no significant difference in BMI among the four groups (Fig 3B bottom). We next visualized how categorical variables differed among the four groups, stratified by gender (Fig 3C). In the figure, 1, 2, 3, 4 refer to groups 1–4 in men, whereas 5, 6, 7, 8 refer to groups 1–4 in women. Group 4 showed higher frequencies of underlying medical conditions (hypertension, diabetes, rheumatism, heart disease, collagen disease, liver disease) and use of medications (steroids, NSAIDs, immunosuppressants) as well as lower frequencies of adverse reactions to vaccinations (local pain, fever, fatigue, headache, joint pain, nausea).

We further performed multiple regression analysis to see whether each of the three features could be explained by the participants’ demographic and health information (S2 Table). The obtained models for the logarithm of the peak and the duration and logarithm of AUC had R2 values of 0.215, 0.267 and 0.217, respectively. The variables that significantly influenced the peak were age, the interval between vaccinations, dyslipidemia, fever over 37.5 degrees, and gender; for duration, they were age, fever, smoking, steroids, immunosuppressants, and dyslipidemia; and for the AUC, they were age, immunosuppressants, steroids, the interval between vaccinations, dyslipidemia, fever over 37.5 degrees, kidney disease, and rheumatism. Hence, all three features were partially explained by the basic demographic and health information. We note that blood groups (A, B, O, AB) and BCG vaccination history were not predictors of these features [27,28]. In fact, there were no significant differences between BCG vaccinated and unvaccinated individuals in peak (p = 0.892), duration (p = 0.521) and area under the curve (p = 0.873). There were also no significant differences between blood groups in peak (p = 0.850), duration (p = 0.211) and area under the curve (p = 0.200).

Deriving a personalized antibody score

Recently, a systematic approach to fit optimized scores with mixed-integer nonlinear programming was proposed [29]. Combining the demographic and health information in Table 1, we devised a simple score aimed at roughly estimating their antibody status. We chose the AUC as a representative feature of individual antibody status because it reflects both the early and the late phases of antibody dynamics. We here constructed two types of scores to cover the whole range of AUC as follows: a score to predict whether an individual’s AUC is in the top third of the population (i.e., top AUC score, Fig 4A) and a score to predict whether an individual’s AUC is in the bottom third (i.e., bottom AUC score, S3 Fig A).

Fig 4. Scoring vaccine-elicited antibody dynamics: (A) Metric for calculating the “top AUC score,” i.e., a score to identify individuals with AUC in the top third of the population, is calculated. (B) Distribution of the top AUC scores in the test dataset is shown: 35, 126, 106, 119, 72, and 8 individuals had scores of -1 or less, 0, 1, 2, 3, and 4, respectively. Those in the top third of the test dataset are shown in yellow, and those not in the top third are shown in orange. The ratio of individuals with AUC in the top third of the test dataset increased as the top AUC score increased. (C) The average AUC tended to increase as the top AUC score increased. (D) The top AUC score was correlated with the IgG(S) titers 90 or 180 days from the second vaccination, calculated from the mathematical model.

We first divided our dataset into training and test datasets. The training dataset consisted of participants in Minami Soma City and Hirata Village (1935 individuals), and the test dataset consisted of participants in Soma City (466 individuals). Six individuals who did not fully answer the questionnaire were excluded from the datasets. We used the training dataset for fitting the scores. The algorithm searches the space of linear combinations of factors with integer coefficients from -5 to +5 to find the best combination to differentiate the population in the top third of the training dataset from the rest (or the bottom third from the rest) (Fig 4A). We then assessed the performance of these scores in the test dataset, that is, we tested whether the scores just created could differentiate the population in the top third of the test dataset from the rest (or the bottom third from the rest).

The top AUC scores in the test dataset were between -1 and 4, except for two individuals with scores of -6 and -5, respectively (Fig 4B and S3 Table). The 8 individuals with scores of 4 (shown in yellow and orange) included 6 individuals (shown in yellow) whose AUC belonged to the top third of the test dataset population (75.0%). The 72 individuals with a score of 3 included 43 individuals with AUC in the top third (59.7%). The 119 individuals with a score of 2 included 50 individuals with AUC in the top third (42.0%). The 106 individuals with a score of 1 included 20 individuals with AUC in the top third (18.9%). The 126 individuals with a score of 0 included 19 individuals with AUC in the top third (15.1%). The 35 individuals with a score of -1 or less included 5 individuals with AUC in the top third (14.3%). Thus, the higher an individual’s top AUC score, the more likely they were to belong to the top third of the population and to have a higher AUC as well (Fig 4C). For example, if a healthy 30-year-old male who drinks two bottles of beer every day had a fever of 37.7 degrees but no nausea after vaccination, his score would be calculated as 1+1 = 2, indicating that his AUC is in the top third with a probability of 42.0%. The top AUC score was also well correlated with IgG(S) titers 90 or 180 days from the second vaccination (Fig 4D). When the threshold is set to 1.5, the top AUC score has an accuracy of 69.1%, a precision of 49.7%, and a recall of 69.2%, evaluated on the test dataset (see Methods for calculation methods).

On the other hand, we calculated the bottom AUC scores in the test dataset, which were between -4 and 2, except for one individual with a score of 3 (S3 Fig B and S3 Table). We confirmed similar trends. For example, the 58 individuals with a score of 1 included 40 individuals with AUC in the bottom third (69.0%). Thus, the higher an individual’s bottom AUC score, the more likely they were to belong to the bottom third of the population and have a lower AUC as well (S3 Fig C). The bottom AUC score was also inversely correlated with IgG(S) titers 90 or 180 days from the second vaccination (S3 Fig D). When the threshold is set to -1.5, the bottom AUC score has an accuracy of 61.2%, a precision of 50.7%, and a recall of 80.5%, evaluated on the test dataset. Note that these AUC scores are useful for detecting the top and bottom populations rather than the middle population, given the relatively low frequency of some of the items listed (collagen disease, immunosuppressant, rheumatism). In fact, when we created a "middle AUC score" to detect individuals with AUC in the middle third (S10 Fig), 92.7% of individuals in the test dataset had the same score of 0 (S3 Table), suggesting that this middle AUC score is not very useful. Note also that even the top and bottom AUC scores will not be able to detect all individuals with top or bottom titer levels, and there will be false positives and false negatives as shown in S3 Table.

We also created a scoring for groups 1 (with high peak and long duration) and 4 (with low peak and short duration) in Fig 3A. Group 1 score had similar items to the top AUC score, and group 4 score had similar items to the bottom AUC score (S11 Fig). These results show that the personalized AUC score can provide a reasonable estimate of an individual’s antibody status, especially for the top and bottom populations, allowing the individual to make informed decisions about disease prevention.

Discussion

In this study, we created a personalized antibody score that can be used as a basis for identifying populations with low sustained antibody titers. To derive an optimal score, we used a mathematical model of antibody production in response to two-dose mRNA vaccinations and reconstructed the vaccine-elicited antibody dynamics of 2,407 participants from the Fukushima vaccination cohort. In particular, we highlighted that the reconstructed best-fit antibody titer curves perfectly predicted the additional timepoint data sampled from 120 of the 2,407 participants. Our mechanism-based mathematical modeling, in contrast to the statistical modeling used in recent reports [30,31], enables a biologically accurate description and precise comparison of antibody dynamics. The parameters of the estimated dynamics showed a large variation spanning two orders of magnitude. This variation was partially explained by individual characteristics like age, sex, the interval between the two vaccine doses, adverse reactions, comorbidities, and medications taken. This result is consistent with previous studies reporting age, sex, vaccine interval, and comorbidities as factors affecting antibody titers [5,32].

Our antibody scores can be easily calculated from individual demographic and health information, yet identify participants with high and low antibody titers (AUC of the IgG(S) titers) to a certain extent. Given the pleiotropic aspect of humoral immunity, it is surprising that a score consisting of several questionnaire items can provide predictions. The score showed that COVID-19 infection history and adverse reactions positively affect the AUC, whereas age and immunologically compromised conditions negatively affect the AUC. These positive and negative factors are consistent with previous studies [5,19,33–35]. In addition, autoimmune diseases have been reported to be risk factors for breakthrough infections [36,37]. Considering that individual antibody titers are partially predictive of the likelihood of breakthrough infections [11,38], this suggests that our antibody score has potential to be used as a risk score for breakthrough infections. On the other hand, the antibody score also has some similarities with COVID-19 severity scores [39–43]. In fact, both include age as a factor likely to lead to low antibody titers and critical illness, both of which may be related to a defective immune system, as observable in cytokine signatures [44,45] or immunoglobulin responses [46,47]. However, whereas COVID-19 severity scores use the results of laboratory tests and clinical symptoms to assess the patient’s condition in the hospital, our antibody score can be easily calculated on the basis of questionnaire responses provided by individuals (i.e., not limited to patients) themselves.

So far, 13.01 billion doses of COVID-19 vaccines have been administered globally, and 2.51 million are now administered each day [48]. Whereas 68.5% of the world population has received at least one dose of a COVID-19 vaccine, 24.6% of people in low-income countries have received at least one dose as the result of disparities in global access to vaccines [49,50]. For example, only a small percentage of persons have received one vaccine dose in African countries [51]. Although the third and fourth doses of COVID-19 vaccines have already begun to be offered, vaccine uptake after the primary two doses is challenged by (booster) vaccine hesitancy, waning of vaccine coverage due to painful adverse reactions [49], and a drop-off in interest in COVID-19 [52]. Even during the unprecedented but successful worldwide vaccination programs for COVID-19, we experienced difficulty in the global distribution of vaccine doses. One of the lessons we learned was the importance of identifying vulnerable populations to achieve early control of infectious disease. In future pandemics, the score that we have developed may play an important role in identifying vulnerable populations.

There are limitations to this study that we can improve for a better understanding of the role of vaccination. First, we measured IgG(S) and neutralization activity against the Wuhan strain as the humoral immune status, and analyzed IgG(S) because of the high correlation between these measurements. Although it is reported that antibodies induced by infection with the ancestral SARS-CoV-2 strain and/or the first of the primary two-dose vaccinations show cross-neutralization of variants from Alpha to Omicron BA.1 [53,54], further analysis of antibody responses (in particular, for neutralization) against variants of concern may be important to extend our scoring to prediction of vaccine efficacy for corresponding strains. Although our current score is limited in its detection capacity, the creation of such a personalized antibody score in the future by using individual-level demographic and health information will be an important tool for designing optimal vaccination strategies in future pandemics (not only COVID-19 but also other infectious diseases). Second, although one of the advantages of the Fukushima vaccination cohort is the extremely low number of natural infections, and that we can evaluate antibody responses minimizing the effect of infection, large portions of the population have been infected with COVID-19 in many parts of the world. In addition, the third and fourth doses of COVID-19 vaccines have already begun to be administered. Our primary purpose in this study was to establish an integrated framework for scoring, but an extended score considering breakthrough infection and booster vaccinations is worth evaluating. Another limitation is that the model fitting was based on antibody measurements and the antibody score was solely based on information available from the questionnaire. This is reflected in the rather low R2 values (around 0.2–0.3) of the regression models of the antibody features. Therefore, further refinement of the score using additional information will be worthwhile. It is worth mentioning that the immune system is affected by multiple factors, including genetics, the environment (such as cohabitation), and markers of metabolic health [55–58], all of which likely influence individual antibody status but were not considered here. The biological determinants of antibody variation will be further revealed in future studies addressing not only B cell subsets but the whole immune system encompassing adaptive as well as innate immunity.

In conclusion, quantifying the variability in antibody dynamics can be a basis for policy decisions regarding the distribution of booster vaccines to strengthen immunity [59] or the use of oral antiviral drugs for the treatment of breakthrough infections [60]. As we learned from the COVID-19 pandemic, determining the priority for (booster) vaccines is an important public health concern given a limited global supply of vaccines at the beginning of a pandemic. For example, since we have data on what percentage of people with certain top/bottom AUC scores have low antibody levels, we can estimate how many people should be prioritized for booster vaccination based on these scores. Thus, the score developed here can be used to estimate and allocate necessary vaccination resources to “priority groups” for vaccination [61,62] in addition to the standard prioritized population [1,2,5,6]. Moreover, the score may be used by medical practitioners to encourage individuals with low predicted antibody levels to get (booster) vaccinations. However, we caution that our current score does not identify all individuals who need booster vaccinations and should not be the sole resource for guiding individual decisions. Taking advantage of our Fukushima vaccination cohort, we will further evaluate the impact of booster vaccinations and post-vaccination infections on personalized antibody scores. It is important that we make use of the lessons learned from COVID-19 vaccination to combat future infectious diseases.

Methods

Ethics statement

The study was approved by the ethics committees of Hirata Central Hospital (number 2021-0611-1) and Fukushima Medical University School of Medicine (number 2021–116). Written informed consent was obtained from all participants individually before the survey.

Participant recruitment and sample collection

The candidates were mainly recruited from Hirata village, Soma city, and Minamisoma city in rural Fukushima prefecture. We conducted non-sequential blood sampling and sequential blood sampling. A total of 2526 individuals participated in non-sequential blood sampling (Fig 1E). Health care workers, frontline workers, and administrative officers from each municipality were intentionally recruited to keep the cohort size large and the dropout rate low. Although most of the health care workers, frontline workers, and administrative officers were under the age of 65, relatively healthy community-dwelling older adults living in the community and in long-term care facilities were also recruited to maintain a wide age range for the cohort. Blood sampling was performed once during each period in June, September, and November 2021, respectively. The first vaccine dose was administered between March 10 and August 20, 2021, and the second dose between March 31 and September 14. The median (interquartile range) interval for the two-dose vaccination was 21 days. A total of 226 health care workers participated in the first blood sampling between May 31 and June 6, 2021. A total of 2526 individuals participated in the second blood sampling between September 8 and October 8, 2021. A total of 2443 individuals participated in the third blood sampling between November 21 and December 25, 2021. In contrast, 12 health care workers participated in the sequential blood sampling, and their vaccination and blood sampling schedule are shown in Fig 1E. Note that the 12 health care workers (described in S5 Fig A) with sequential blood sampling were not included in the non-sequential population of 2526 participants. In conclusion, of the total 2526 participants, those eligible for analysis were 2459 participants who completed the second vaccination and at least two blood samplings.

Out of these 2459 participants, 43 individuals had higher IgG(S) titers on the second measurement than on the first, although their IgG(N) titers were negative (meaning no infection history). We speculated that this was due to measurement error and did not attempt model fitting. Furthermore, 9 individuals who had an interval of longer than 40 days between the first and the second vaccination were also removed. Because most of our participants had an interval of 21 or 28 days and there were very few cases of longer intervals, we decided that there were not enough data to reliably perform predictions among the participants with longer intervals. This was a conservative decision to ensure the high quality of our model fitting. As a result, we performed parameter estimation on the remaining 2407 participants (Fig 1B).

Information on sex, age, daily medication, medical history, date of vaccination, adverse reaction after vaccination, type of vaccination, blood type, bacillus Calmette–Guérin (BCG) vaccine history, smoking habits, and drinking habits was retrieved from the paper-based questionnaire (summarized in Table 1). All blood sampling was performed at the medical facilities with 8 mL, and serum samples were sent to The University of Tokyo.

SARS-CoV-2-specific antibody measurement

All serological assays were conducted at The University of Tokyo. Specific IgG (i.e., IgG(S)) and neutralizing activity against the Wuhan strain were measured as the humoral immune status after COVID-19 vaccination. Specific IgG antibody titers (IgG(N)) were used to determine past COVID-19 infection status. Chemiluminescent immunoassay with iFlash 3000 (YHLO Biotech, Shenzhen, China) and iFlash-2019-nCoV series (YHLO Biotech, Shenzhen, China) reagents were used in the present study. The measurement range was 2–3500 AU/mL for IgG(S) and 4–800 AU/mL for neutralizing activity. For neutralizing activity, AU/mL×2.4 was used to convert to International Units (IU/mL); for IgG(S), AU/mL×1.0 was used to convert to binding antibody units (BAU/mL). The testing process was as per the official guideline. Quality checks were conducted every day before starting the measurement.

Modeling vaccine-elicited antibody dynamics

We developed a mathematical model describing COVID-19 vaccine-elicited antibody dynamics to evaluate the impact of primary two-dose COVID-19 vaccination on rapid immunity at the individual level and reconstructed the best-fit antibody titer curves of 2,407 participants in the Fukushima vaccination cohort. Here we explain the derivation and formulation of the mathematical model in detail.

(i) Vaccination-elicited antibody dynamics after the first dose.

After the first vaccination, naïve B cells encounter antigens and differentiate into short-lived antibody-secreting cells (ASCs), plasmablasts, germinal center (GC) B cells, or GC-independent memory B cells depending on BCR affinity for their cognate antigen [63]. Then, the GC B cells undergo rapid proliferation with somatic immunoglobulin hypermutation and subsequently differentiate into GC-dependent memory B cells or long-lived antibody-secreting cells, which are plasma cells with immunoglobulin class switching. To describe this antigen-specific B cell expansion and the induction of antibody-secreting cells and memory B cells after the first vaccination (S4 Fig A), we developed a simple but quantitative mathematical model as follows:[Equation Omitted: See PDF][Equation Omitted: See PDF][Equation Omitted: See PDF]where the variables M(t) = M1(t), B(t), and A(t) are the amount of mRNA inoculated by the vaccination, the number of antibody-secreting cells, and the antibody titers at time t, respectively. The parameters τ1 and d represent the timing of the vaccination and the decay rate of mRNA. We here considered D1 to be the inoculated dose of mRNA by the vaccination, that is, M1(τ1) = D1.

Because the data we used here were limited (i.e., only time-course vaccine-elicited IgG(S) titers), one compartment of B cells including heterogeneous cell populations that produce antibodies (i.e., short-lived and long-lived antibody-secreting cells) was assumed. Therefore, we modelled the average B cell population dynamics in Eq (2), where the product of P1(t) and M(t)m/(M(t)m + Km) represents the average de novo induction of the antibody-secreting cells. P1(t) is a step function defined as P1(t) = P1 for τ1 + η1 ≤ t, where η1 is the delay of induction of antibody-secreting cells after vaccination: otherwise P1(t) = 0. The parameters m, K, and μ correspond to the steepness at which the induction increases with increasing amount of mRNA (i.e., the Hill coefficient), the amount of mRNA satisfying P1/2, and the average decay rate of the antibody-secreting cell compartment, respectively. The other parameters, p and c, represent the antibody production rate and the clearance rate of antibodies, respectively.

(ii) Vaccination-elicited antibody dynamics after the second dose.

After the second vaccination, memory B cells are reactivated by re-exposure to the antigen. Some differentiate into short-lived antibody-secreting cells (plasmablasts) or memory B cells outside the GC. Others enter the GC to become secondary GC B cells. Subsequently, these secondary GC B cells differentiate into GC-dependent memory B cells or long-lived antibody-secreting cells (plasma cells). To describe these recall B cell responses and their secretion of antibody after the second vaccination (S4 Fig B), we modified the above mathematical model, Eqs (1–3), as follows:[Equation Omitted: See PDF][Equation Omitted: See PDF][Equation Omitted: See PDF]where M2(t) is the amount of mRNA by the vaccination inoculated at τ2 satisfying M2(τ2) = D2. In addition, P2(t) = P1 and M(t) = M1(t) for τ2≤t<τ2+η2, P2(t) = P2 and M(t) = M1(t) + M2(t) for τ2 + η2 ≤ t, where η2 is the delay of induction of antibody-secreting cells after vaccination; otherwise P2(t) = 0. Here we ignored the GC-independent memory B cells induced by the first vaccination because of their minor effect. In general, once reactivated, memory B cells can re-enter the GC more rapidly than naïve B cells, and therefore the secondary antibody responses are much faster and larger (i.e., η1 >η2 and P1 < P2, respectively). In the main recall immunity, the quantity and quality of memory B cells established by the first vaccination is included in P2.

(iii) Mathematical model for data fitting.

Since the clearance rate of antibody is much larger than the decay of antibody-secreting cells (i.e., c ≫ μ), we made a quasi-steady state assumption, dA(t)/dt = 0 dV(t)/dt = 0, and replaced Eqs (3) and (6) with A(t) = pB(t)/c. Moreover, since Eqs (1) and (4) are linear differential equations, M1(t) = D1e-dt for t ≥ τ1 and M2(t) = D2e-dt for t ≥ τ2: otherwise M1(t) = M2(t) = 0, respectively. Thus, the above Eqs (1–6) are further simplified assuming τ1 = 0 and D1 > 0, and we obtained the following single ordinary differential equation, which we used to analyze the antibody responses (i.e., IgG(S) titers (AU/mL)) in this study:[Equation Omitted: See PDF]Where H(t) = Hi = pPi/c for τi + ηi ≤ t<τi+1 + ηi+1 (i = 1 or 2), η3 = ∞ and D2 > 0 for τ2 + η2 ≤ t: otherwise D2 = 0. This simple model can quantify the vaccine-elicited time-course antibody dynamics as described in Fig 2A under an arbitrary threshold of antibody titers ATH (see below).

Quantifying vaccine-elicited time-course antibody dynamics

In addition to the participants in the cohort, we included 12 health care workers whose serum was sequentially sampled for 40 days (on average 25 samples per individual) for validation and parameterization of a mathematical model for vaccine-elicited antibody dynamics. A nonlinear mixed effects model was used to fit the antibody dynamics model, given by Eq (7), to the longitudinal antibody titers of IgG(S) obtained from the 12 health care workers. The mathematical model included both a fixed effect and a random effect in each parameter. That is, the parameters for individual k, θk(= θ×eπk) are represented as a product of θ (a fixed effect) and eπk (a random effect). πk follows a normal distribution with mean 0 and standard deviation Ω. As we described above, since η1 > η2 and P1 < P2, we assumed η1 = η and η2 = fdelayη, H2 = H and H1 = fdegreeH, and estimated η, H, 0 < fdelay < 1 and 0 < fdegree < 1 for conducting biologically reasonable estimations. We here assumed that the parameters η, H, fdelay, fdegree and m varied across individuals, whereas we did not consider interindividual variability in other parameters to ensure parameter identifiability. Note that the half-life of mRNA (i.e., log 2/d) and dose of mRNA (i.e., Di) are assumed to be 1 day [64] and 100 (μg/0.5 mL) [65], respectively. Fixed effect and random effect were estimated by using the stochastic approximation expectation-approximation algorithm and empirical Bayes’ method, respectively. Fitting was performed using MONOLIX 2019R2 (www.lixoft.com) [66]. The estimated (fixed and individual) parameters are listed in S1 Table. With the estimated parameters for each individual, the dynamics of IgG(S) titers, A(t), and the average de novo antibody response elicited by the first and second vaccinations, , were calculated in S5 Fig AB, respectively. Interestingly, we observed that the variations induced by the second vaccination were much larger than those induced by the first vaccination.

We found that most of the best-fitted estimated parameters in the mathematical model (i.e., D1, D2, d, μ, K, η, fdelay, fdegree) were the same or similar across the 12 individuals compared with those of parameters of m and H (see S1 Table). We note that m and H, which showed wide variation of estimated values, contributed mainly to the vaccine-elicited antibody dynamics after the second vaccine dose, whereas the other parameters contributed to that after the first dose. In fact, we are interested in the large variation after the second dose (but not the negligible variation after the first dose). Therefore, we hereafter fixed the parameters in our mathematical model to be the estimated population parameters listed in S1 Table, except m and H. These assumptions enabled us to accurately reconstruct the large variations in antibody dynamics after the second dose, and these two parameters were independently estimated from each IgG(S) by a nonlinear least-squares method, even if 2,407 participants had only 2 or 3 measurements of antibody titers at different time points. Using the estimated parameters for each individual, we fully reconstructed the dynamics of IgG(S) titers of all participants after the first vaccination in S5 Fig C. We summarized the distribution of parameter values m and H for 2,407 participants in S5 Fig D, and the best-fit antibody titer curves of 200 randomly selected individuals are plotted along with the observed data for visualization in S2 Fig.

Mathematical model validation

Outside the study period described in Fig 1E, to validate our mathematical model (i.e., Eqs (4–5)), we first prepared Validation dataset A: the additional sequentially performed serum sampling of 12 health care workers, that is, data for 3 to 4 timepoints before the booster dose (on average 4 additional samples per individual for 165 days after the first vaccine dose). Employing our estimated parameters, which were estimated from the original dataset (i.e., the black circles in S6 Fig), listed in S1 Table, we found our mathematical model predicted the values for the additional 3 to 4 timepoints before booster vaccination with high accuracy. Next, to validate the reconstructed best-fit antibody titer curves described in S5 Fig C, we also prepared Validation dataset B: an additional third or fourth antibody measurement from 110 of the 2407 participants of our Fukushima vaccination cohort around 3 to 4 months after the second or third timepoint data but before the booster vaccination. We compared the reconstructed antibody titer curves and the additional data and confirmed that the additional data points obtained from participants who were not infected with COVID-19 (i.e., IgG(N)-negative participants) during that period perfectly matched the prediction of our mathematical model (S7 Fig). On the other hand, interestingly, for those who got infected with COVID-19, the additional data were significantly larger than the reconstructed curves (S8 Fig).

Additionally, we assessed the accuracy of the reconstructed best-fit antibody titers by calculating the Pearson’s correlation coefficient between the reconstructed antibody titers and the observed (or additionally observed) antibody titers (S9 Fig). The Pearson’s correlation coefficients between the observed antibody titers, which were used in parameter estimation, and the reconstructed antibody titers, were 0.99 for both of the 12 health care workers and the 2407 participants in the Fukushima vaccination cohort. Furthermore, the Pearson’s correlation coefficients between the additionally observed antibody titers, which were not used in parameter estimation, and the reconstructed antibody tiers, were 0.94 and 0.92 for the 12 health care workers and the 110 participants in the Fukushima vaccination cohort, respectively.

Taken together, our additional data (i.e., Validation dataset A and B) and quantitative analysis demonstrated that our mathematical model and the reconstructed best-fit antibody titer curves using the same are validated for the purposes of predicting antibody dynamics in a generalized population.

Building optimized antibody scores

A Python implementation of Ustun et al.’s [29] algorithm (risk-slim, https://github.com/ustunb/risk-slim) was used to build optimized AUC scores. Briefly, the algorithm searches for the best linear combination of features with integer coefficients that minimizes the sum of the logistic loss and the l0-norm of the coefficients. The range of coefficients was set to -5 to +5; the l0-penalty parameter C0 was set to 1×10−3.

Evaluation of antibody scores

When an individual’s top (or bottom) AUC score is above a threshold, that individual is predicted to be in the top 1/3 (or the bottom 1/3) of the population. The accuracy, precision, and recall of the score are defined as (TP + TN)/(TP + TN + FP + FN), TP/(TP + FP), and TP/(TP + FN), respectively, where TP, TN, FP, and FN denote true positives (i.e., the number of individuals predicted to be positive that were actually positive), true negatives, false positives, and false negatives, respectively.

Statistical analysis

The answers from the 2,407 participants who completed the paper-based questionnaire were converted into a set of categorical and numerical variables. Numerical variables included age, BMI, and the interval between the two doses. These variables were then used in a multiple regression analysis to explain the three antibody dynamics features. Six participants who did not fully answer the questionnaire were excluded from the analysis. The variables used here belonged to any of the five categories: (i) basic demographic information and lifestyle habits, (ii) information on vaccinations, (iii) underlying medical conditions, (iv) adverse reactions, and (v) medications being taken. When necessary, the same variables were compared among different generations or different groups using Pearson’s chi-square test (for categorical variables), analysis of variance (ANOVA, for more than two numerical variables), or Welch T-test (for two numerical variables). A Bonferroni correction was applied for multiple comparisons. All statistical analyses were performed using R (version 4.1.2) or JMP Pro 16.

Map of Japan

The basemap shapefile of Japan and Fukushima Prefecture was downloaded from National Land Numerical Information (https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_4.html) and edited using the R package sf.

Supporting information

S1 Fig. Comparison between IgG(S) and neutralization activity.

https://doi.org/10.1371/journal.pdig.0000497.s001

(DOCX)

S2 Fig. Reconstructed antibody titer trajectory for individual participants.

https://doi.org/10.1371/journal.pdig.0000497.s002

(DOCX)

S3 Fig. Analyzing antibody titers.

https://doi.org/10.1371/journal.pdig.0000497.s003

(DOCX)

S4 Fig. Modeling vaccine-elicited B cell dynamics.

https://doi.org/10.1371/journal.pdig.0000497.s004

(DOCX)

S5 Fig. Calibrating vaccine-elicited antibody dynamics.

https://doi.org/10.1371/journal.pdig.0000497.s005

(DOCX)

S6 Fig. Validating calculated lgG(S) titers for each individual health care worker.

https://doi.org/10.1371/journal.pdig.0000497.s006

(DOCX)

S7 Fig. Validating reconstructed lgG(S) titers for 109 participants from the Fukushima vaccination cohorts who did not receive the booster vaccination and did not get infected with COVID-19.

https://doi.org/10.1371/journal.pdig.0000497.s007

(DOCX)

S8 Fig. Validating reconstructed lgG(S) titers for 10 participants from the Fukushima vaccination cohorts who did not receive the booster vaccination and got infected with COVID-19.

https://doi.org/10.1371/journal.pdig.0000497.s008

(DOCX)

S9 Fig. Evaluating estimation accuracy by mathematical model through a comparison of observed and estimated antibody titers.

https://doi.org/10.1371/journal.pdig.0000497.s009

(DOCX)

S10 Fig. Middle AUC scores.

https://doi.org/10.1371/journal.pdig.0000497.s010

(DOCX)

S11 Fig. Scoring for groups 1 and 4.

https://doi.org/10.1371/journal.pdig.0000497.s011

(DOCX)

S1 Table. Estimated fixed and individual parameters for 12 health care workers.

https://doi.org/10.1371/journal.pdig.0000497.s012

(DOCX)

S2 Table. Coefficients and p-values of the multiple regression analysis.

https://doi.org/10.1371/journal.pdig.0000497.s013

(DOCX)

S3 Table. Distributions of the top, bottom and middle AUC scores.

https://doi.org/10.1371/journal.pdig.0000497.s014

(DOCX)

S1 Data. Synthetic patient dataset.

https://doi.org/10.1371/journal.pdig.0000497.s015

(CSV)

Acknowledgments

We would like to thank all the staff from Fukushima Medical University, Seireikai health care group, Hirata Village office, Soma City office, Soma Central Hospital, Soma General Hospital, Minamisoma City office, Minamisoma City Medical Association, Minamisoma Municipal General Hospital, Shindo Clinic and Medical Governance Institute, who contributed significantly to the accomplishment of this research, especially, Ms. Yuka Harada, Ms. Serina Noji, Ms. Naomi Ito, Dr. Makoto Kosaka, Mr. Anju Murayama, Mr. Sota Sugiura, Mr. Manato Tanaka, Ms. Yuna Uchi, Mr. Yudai Kaneda, Mr. Masahiko Nihei, Mr. Hideo Sato, Ms. Rie Yanai, Ms. Yasuko Suzuki, Ms. Keiko Abe, Dr. Hidekiyo Tachiya, Mr. Kouki Nakatsuka, Dr. Ryuzaburo Shineha, Ms. Miki Sato, Dr. Masahiko Sato, Mr. Naoharu Tadano, Mr. Kazuo Momma, Mr. Shu-ichi Mori, Ms. Saori Yoshisato, Ms. Katsuko Onoda, Mr. Satoshi Kowata, Mr. Masatsugu Tanaki, Dr. Tomoyoshi Oikawa, Dr. Joji Shindo, Ms. Xujin Zhu, Ms. Asaka Saito, Ms. Yuumi Kondo, and Ms. Tomoyo Nishimura.

Citation: Nakamura N, Kobashi Y, Kim KS, Park H, Tani Y, Shimazu Y, et al. (2024) Modeling and predicting individual variation in COVID-19 vaccine-elicited antibody response in the general population. PLOS Digit Health 3(5): e0000497. https://doi.org/10.1371/journal.pdig.0000497

References

1. Goldberg Y, Mandel M, Bar-On YM, Bodenheimer O, Freedman L, Haas EJ et al., Waning Immunity after the BNT162b2 Vaccine in Israel. N Engl J Med 385, e85 (2021). pmid:34706170

2. Levin EG, Lustig Y, Cohen C, Fluss R, Indenbaum V, Amit S et al., Waning Immune Humoral Response to BNT162b2 Covid-19 Vaccine over 6 Months. N Engl J Med 385, e84 (2021). pmid:34614326

3. Shroff RT, Chalasani P, Wei R, Pennington D, Quirk G, Schoenle MV et al., Immune responses to two and three doses of the BNT162b2 mRNA vaccine in adults with solid tumors. Nat Med 27, 2002–2011 (2021). pmid:34594036

4. Choi A, Koch W, Wu K, Chu L, Ma LZ, Hill A et al., Safety and immunogenicity of SARS-CoV-2 variant mRNA vaccine boosters in healthy adults: an interim analysis. Nat Med 27, 2025–2031 (2021). pmid:34526698

5. Ward H, Whitaker M, Flower B, Tang SN, Atchison C, Darzi A et al., Population antibody responses following COVID-19 vaccination in 212,102 individuals. Nat Commun 13, 907 (2022). pmid:35173150

6. Lustig Y, Sapir E, Regev-Yochay G, Cohen C, Fluss R, Olmer L et al., BNT162b2 COVID-19 vaccine and correlates of humoral immune responses and dynamics: a prospective, single-centre, longitudinal cohort study in health-care workers. Lancet Respir Med 9, 999–1009 (2021). pmid:34224675

7. Yau K, Abe KT, Naimark D, Oliver MJ, Perl J, Leis JA et al., Evaluation of the SARS-CoV-2 Antibody Response to the BNT162b2 Vaccine in Patients Undergoing Hemodialysis. JAMA Netw Open 4, e2123622 (2021). pmid:34473256

8. Moyon Q, Sterlin D, Miyara M, Anna F, Mathian A, Lhote R et al., BNT162b2 vaccine-induced humoral and cellular responses against SARS-CoV-2 variants in systemic lupus erythematosus. Ann Rheum Dis 81, 575–583 (2022). pmid:34607791

9. Lumley SF, O’Donnell D, Stoesser NE, Matthews PC, Howarth A, Hatch SB et al., Antibody Status and Incidence of SARS-CoV-2 Infection in Health Care Workers. N Engl J Med 384, 533–540 (2021). pmid:33369366

10. Williams DE, Statistics of antibody binding to the spike protein explain the dependence of COVID 19 infection risk on antibody concentration and affinity. Sci Rep 12, 9379 (2022). pmid:35672372

11. Lipsitch M, Krammer F, Regev-Yochay G, Lustig Y, Balicer RD, SARS-CoV-2 breakthrough infections in vaccinated individuals: measurement, causes and impact. Nat Rev Immunol 22, 57–65 (2022). pmid:34876702

12. Kobashi Y, Kawamura T, Shimazu Y, Zhao T, Sugiyama A, Nakayama A et al., Humoral immunity after second dose of BNT162b2 vaccine in Japanese communities: an observational cross-sectional study, Fukushima Vaccination Community Survey. Sci Rep 12, 18929 (2022). pmid:36344597

13. Kobashi Y, Nishikawa Y, Kawamura T, Kodama T, Shimazu Y, Obara D, Zhao T, Tsubokura M et al., Seroprevalence of SARS-CoV-2 antibodies among hospital staff in rural Central Fukushima, Japan: A historical cohort study. Int Immunopharmacol 98, 107884 (2021). pmid:34246041

14. Kobashi Y, Shimazu Y, Kawamura T, Nishikawa Y, Omata F, Kaneko Y et al., Peak IgG antibody titers against SARS-CoV-2 spike protein following immunization with the Pfizer/BioNTech BNT162b2 vaccine. Fukushima J Med Sci 68, 67–70 (2022). pmid:35228456

15. Kobashi Y, Shimazu Y, Kawamura T, Nishikawa Y, Omata F, Kaneko Y et al., Factors associated with anti-severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) spike protein antibody titer and neutralizing activity among healthcare workers following vaccination with the BNT162b2 vaccine. PLoS One 17, e0269917 (2022). pmid:35687563

16. Yoshida M, Kobashi Y, Kawamura T, Shimazu Y, Nishikawa Y, Omata F et al., Factors Associated with COVID-19 Vaccine Booster Hesitancy: A Retrospective Cohort Study, Fukushima Vaccination Community Survey. Vaccines (Basel) 10 (2022). pmid:35455264

17. Yoshida M, Kobashi Y, Shimazu Y, Saito H, Yamamoto C, Kawamura T et al., Time course of adverse reactions following BNT162b2 vaccination in healthy and allergic disease individuals aged 5–11 years and comparison with individuals aged 12–15 years: an observational and historical cohort study. Eur J Pediatr , 1–11 (2022). pmid:36224435

18. Zhao T, Nishi-Uchi T, Omata F, Takita M, Kawashima M, Nishikawa Y et al., Humoral response to SARS-CoV-2 vaccination in haemodialysis patients and a matched cohort. BMJ Open 12, e065741 (2022). pmid:36351730

19. Tani N, Chong Y, Kurata Y, Gondo K, Oishi R, Goto T et al., Relation of fever intensity and antipyretic use with specific antibody response after two doses of the BNT162b2 mRNA vaccine. Vaccine 40, 2062–2067 (2022). pmid:35177298

20. Hall VJ, Foulkes S, Saei A, Andrews N, Oguti B, Charlett A et al., COVID-19 vaccine coverage in health-care workers in England and effectiveness of BNT162b2 mRNA vaccine against infection (SIREN): a prospective, multicentre, cohort study. Lancet 397, 1725–1735 (2021). pmid:33901423

21. Rhodes S, Wilkinson J, Pearce N, Mueller W, Cherrie M, Stocking K et al., Occupational differences in SARS-CoV-2 infection: analysis of the UK ONS COVID-19 infection survey. J Epidemiol Community Health 76, 841–846 (2022). pmid:35817467

22. Shrotri M, Krutikov M, Palmer T, Giddings R, Azmi B, Subbarao S et al., Vaccine effectiveness of the first dose of ChAdOx1 nCoV-19 and BNT162b2 against SARS-CoV-2 infection in residents of long-term care facilities in England (VIVALDI): a prospective cohort study. Lancet Infect Dis 21, 1529–1538 (2021). pmid:34174193

23. Xie J, Feng S, Li X, Gea-Mallorquí E, Prats-Uribe A, Prieto-Alhambra D et al., Comparative effectiveness of the BNT162b2 and ChAdOx1 vaccines against Covid-19 in people over 50. Nat Commun 13, 1519 (2022). pmid:35314696

24. Cromer D, Steain M, Reynaldi A, Schulub TE, Wheatley AK, Juno JA et al., Neutralising antibody titres as predictors of protection against SARS-CoV-2 variants and the impact of boosting: a meta-analysis. Lancet Microbe 3, e52–e61 (2022). pmid:34806056

25. Khoury DS, Cromer D, Reynaldi A, Schulub TE, Wheatley AK, Juno JA et al., Neutralizing antibody levels are highly predictive of immune protection from symptomatic SARS-CoV-2 infection. Nat Med 27, 1205–1211 (2021). pmid:34002089

26. Goel RR, Painter MM, Lundreen KA, Apostolidis SA, Baxter AE, Giles JR et al., Efficient recall of Omicron-reactive B cell memory after a third dose of SARS-CoV-2 mRNA vaccine. Cell 185, 1875–1887 (2022). pmid:35523182

27. Gonzalez-Perez M, Sanchez-Tarjuelo R, Shor B, Nistal-Villan E, Ochando J, The BCG Vaccine for COVID-19: First Verdict and Future Directions. Front Immunol 12, 632478 (2021). pmid:33763077

28. Kander T, Bjurstrom MF, Frigyesi A, Joud M, Nilsson CU, ABO and RhD blood group are not associated with mortality and morbidity in critically ill patients; a multicentre observational study of 29 512 patients. BMC Anesthesiol 22, 91 (2022). pmid:35366803

29. Ustun B, Rudin C, Learning Optimized Risk Scores. J Mach Learn Res 20 (2019).

30. Perez-Alos L et al., Modeling of waning immunity after SARS-CoV-2 vaccination and influencing factors. Nat Commun 13, 1614 (2022). pmid:35347129

31. Wei J, Pouweis KB, Stoesser N, Matthews PC, Diamond I, Studley R et al., Antibody responses and correlates of protection in the general population after two doses of the ChAdOx1 or BNT162b2 vaccines. Nat Med 28, 1072–1082 (2022). pmid:35165453

32. Parry H, Bruton R, Stephens C, Benteley C, Brown K, Amirthalingam G et al., Extended interval BNT162b2 vaccination enhances peak antibody generation. NPJ Vaccines 7, 14 (2022). pmid:35087066

33. Soetedjo NNM, Iryaningrum MR, Lawrensia S, Permana H, Antibody response following SARS-CoV-2 vaccination among patients with type 2 diabetes mellitus: A systematic review. Diabetes Metab Syndr 16, 102406 (2022). pmid:35104750

34. Grainger R, Kim AHJ, Conway R, Yazdany J, Robinson PC, COVID-19 in people with rheumatic diseases: risks, outcomes, treatment considerations. Nat Rev Rheumatol 18, 191–204 (2022). pmid:35217850

35. Kanizsai A, Molnar T, Varnai R, Zavori L, Tőkés-Füzesi M, Szalai Z et al., Fever after Vaccination against SARS-CoV-2 with mRNA-Based Vaccine Associated with Higher Antibody Levels during 6 Months Follow-Up. Vaccines (Basel) 10 (2022). pmid:35335080

36. Ahmed S, Mehta P, Paul A, Anu S, Cherian S, Shenoy V et al., Postvaccination antibody titres predict protection against COVID-19 in patients with autoimmune diseases: survival analysis in a prospective cohort. Ann Rheum Dis 81, 868–874 (2022). pmid:35135831

37. Sun J, Zheng Q, Madhira V, Olex AL, Anzalone AJ, Vinson A et al., Association Between Immune Dysfunction and COVID-19 Breakthrough Infection After SARS-CoV-2 Vaccination in the US. JAMA Intern Med 182, 153–162 (2022). pmid:34962505

38. Tay MZ, Rouers A, Fong S-W, Goh YS, Chan Y-H, Chang ZW et al., Decreased memory B cell frequencies in COVID-19 delta variant vaccine breakthrough infection. EMBO Mol Med 14, e15227 (2022). pmid:34994081

39. Goodacre S, Thomas B, Sutton L, Burnsall M, Lee E, Bradburn M et al., Derivation and validation of a clinical severity score for acutely ill adults with suspected COVID-19: The PRIEST observational cohort study. PLoS One 16, e0245840 (2021). pmid:33481930

40. Liang W, Liang H, Ou L, Chen B, Chen A, Li C et al., Development and Validation of a Clinical Risk Score to Predict the Occurrence of Critical Illness in Hospitalized Patients With COVID-19. JAMA Intern Med 180, 1081–1089 (2020). pmid:32396163

41. Gude F, Riveiro V, Rodríguez-Núñez N, Ricoy J, Lado-Baleato Ó, Lourido T et al., Development and validation of a clinical score to estimate progression to severe or critical state in COVID-19 pneumonia hospitalized patients. Sci Rep 10, 19794 (2020). pmid:33188225

42. Fumagalli C, Rozzini R, Vannini M, Coccia F, Cesaroni G, Mazzeo F et al., Clinical risk score to predict in-hospital mortality in COVID-19 patients: a retrospective cohort study. BMJ Open 10, e040729 (2020). pmid:32978207

43. Huespe I, Bisso IC, Stefano SD, Terrasa S, Gemelli NA, Heras ML, COVID-19 Severity Index: A predictive score for hospitalized patients. Med Intensiva (Engl Ed) (2020). pmid:33478781

44. Del Valle DM, Kim-Schulze S, Huang H-H, Beckman ND, Nirenberg S, Wang B et al., An inflammatory cytokine signature predicts COVID-19 severity and survival. Nat Med 26, 1636–1643 (2020). pmid:32839624

45. Lucas C, Wong P, Klein J, Castro TBR, Silva J, Sundaram M et al., Longitudinal analyses reveal immunological misfiring in severe COVID-19. Nature 584, 463–469 (2020). pmid:32717743

46. Zohar T, Loos C, Fischinger S, Atyeo C, Wang C, Slein MD et al., Compromised Humoral Functional Evolution Tracks with SARS-CoV-2 Mortality. Cell 183, 1508–1519 e1512 (2020). pmid:33207184

47. Lucas C, Klein J, Sundaram ME, Liu F, Wong P, Silva J et al., Delayed production of neutralizing antibodies correlates with fatal COVID-19. Nat Med 27, 1178–1186 (2021). pmid:33953384

48. Our World In Data (2022) Coronavirus (COVID-19) vaccinations. pp https://ourworldindata.org/covid-vaccinations.

49. Arce JSS, Warren SS, Meriggi NF, Scacco A, McMurry N, Voors M et al., COVID-19 vaccine acceptance and hesitancy in low- and middle-income countries. Nat Med 27, 1385–1394 (2021). pmid:34272499

50. VanWormer JJ, Alicea G, Weichelt BP, Berg RL, Sundaram ME, COVID-19 vaccine coverage disparities in rural and farm children. Vaccine (2022). pmid:36400661

51. Lawal L, bello MA, Murwira T, Avoka C, Ma’aruf SY, Omonhinmin IH et al., Low coverage of COVID-19 vaccines in Africa: current evidence and the way forward. Hum Vaccin Immunother 18, 2034457 (2022). pmid:35240908

52. Buneviciene I, Bunevicius R, Bagdonas S, Bunevicius A, COVID-19 media fatigue: predictors of decreasing interest and avoidance of COVID-19-related news. Public Health 196, 124–128 (2021). pmid:34192604

53. Windsor IW, Tong P, Lavidor O, Moghaddam AS, McKay LGA, Gautam A et al., Antibodies induced by an ancestral SARS-CoV-2 strain that cross-neutralize variants from Alpha to Omicron BA.1. Sci Immunol 7, eabo3425 (2022). pmid:35536154

54. Wall EC, Wu M, Harvey R, Kelly G, Warchal S, Sawyer C et al., Neutralising antibody activity against SARS-CoV-2 VOCs B.1.617.2 and B.1.351 by BNT162b2 vaccination. Lancet 397, 2331–2333 (2021). pmid:34090624

55. Liston A, Carr EJ, Linterman MA, Shaping Variation in the Human Immune System. Trends Immunol 37, 637–646 (2016). pmid:27692231

56. Lakshmikanth T, Muhammad SA, Olin A, Chen Y, Mikes J, Fagerberg L et al., Human Immune System Variation during 1 Year. Cell Rep 32, 107923 (2020). pmid:32697987

57. Brodin P, Jojic V, Gao T, Bhattacharya S, Angel CJL, Furman D et al., Variation in the human immune system is largely driven by non-heritable influences. Cell 160, 37–47 (2015). pmid:25594173

58. Liston A, Humblet-Baron S, Duffy D, Goris A, Human immune diversity: from evolution to modernity. Nat Immunol 22, 1479–1489 (2021). pmid:34795445

59. Patel MK, Booster Doses and Prioritizing Lives Saved. N Engl J Med 385, 2476–2477 (2021). pmid:34879192

60. Dal-Re R, Becker SL, Bottieau E, Holm S, Availability of oral antivirals against SARS-CoV-2 infection and the requirement for an ethical prescribing approach. Lancet Infect Dis (2022). pmid:35364019

61. Aukes L, Fireman B, Lewis E, Timbol J, Hansen J, Yu H et al., A Risk Score to Predict Clostridioides difficile Infection. Open Forum Infect Dis 8, ofab052 (2021). pmid:33738316

62. Meier MA, Berger CT, A simple clinical score to identify likely hepatitis B vaccination non-responders ‐ data from a retrospective single center study. BMC Infect Dis 20, 891 (2020). pmid:33238923

63. Laidlaw BJ, Ellebedy AH, The germinal centre B cell response to SARS-CoV-2. Nat Rev Immunol 22, 7–18 (2022). pmid:34873279

64. Anand P, Stahel VP, Review the safety of Covid-19 mRNA vaccines: a review. Patient Saf Surg 15, 20 (2021). pmid:33933145

65. Lin D-Y, Gu Y, Wheeler B, Young H, Holloway S, Sunny S-K et al., Effectiveness of Covid-19 Vaccines over a 9-Month Period in North Carolina. N Engl J Med 386, 933–941 (2022). pmid:35020982

66. Traynard P, Ayral G, Twarogowska M, Chauvin J, Efficient Pharmacokinetic Modeling Workflow With the MonolixSuite: A Case Study of Remifentanil. CPT Pharmacometrics Syst Pharmacol 9, 198–210 (2020). pmid:32036625

67. Japan Geographic Data Center (2021) April 2021 Survey Population and Households by Municipality (Japanese Residents + Foreign Residents) (Japanese). (Japan Geographic Data Center), p https://www.kokudo.or.jp/service/data/map/fukushima.pdf.

68. Fukushima Prefecture (2021) The incidence of COVID-19 in Fukushima Prefecture (Japanese). (Fukushima Prefecture), pp https://www.pref.fukushima.lg.jp/sec/21045c/fukushima-hasseijyoukyou0412.html.