How long someone has detectable neutralizing antibodies after vaccination is important in understanding the impact of vaccination on disease transmission. In this post I step through several different models and re-examine data from a prior paper.

A prominant scientist has been posting his highlights of pre-prints on twitter during the course of the pandemic and occasionally, I come across one that makes me scratch my head. I think that the move to open and rapid science through pre-prints has been good, but not without issues. Review (and open review especially) makes work better and is an important part of science and progress.

One paper that stood out to me was by Suthar et al. (2021) in which the authors look at the neutralizing antibody tires (nAbs) amongst T-cell response in those persons who were vaccinated with the Pfizer SARS-CoV-2 Vaccines (BNT162b2 vaccine). The authors found: “substantial waning of antibody response and T-cell Immunity” in those who were vaccinated. Naturally, this *could* be concerning because antibodies and T-cells are important in fighting infection and could help us understand the endgame of this first epidemic of SARS-CoV-2 (e.g., vaccine induced reduction of transmission, symptomatic disease, severe outcomes, and associated periodicity given time since last vaccine). However, there is the open question of what is the best correlate of protection (e.g., we know that antibodies wane with time, but the ability to produce antibodies in the face of a new infection from memory B-cells is also very important) (Corbett et al. 2021)

As mentioned, along with the bad things about pre-prints (media jumping on sensationalist findings, people seeking acclaim by making bold statements that are not rigorously founded in fact, general review help, etc) is that much of the data are available for secondary analysis. With the article in question, I could find some of the original data analyzed by the authors from a prior publication. We can start to look at these data ourselves:

```
library(data.table)
library(tidyverse)
theme_set(ggdist::theme_ggdist())
dat <- readxl::read_excel("nAbs.xlsx")
dat %>%
gather(day, value, -PID, -Gender) %>%
mutate(day = as.numeric(str_extract(string = day, "\\d+")))->dat_long
knitr::kable(dat, digits = 1)
```

PID | Gender | D_0 | D_21 | D_42 | D_90_120 |
---|---|---|---|---|---|

2000 | Male | 10.0 | 62.3 | 1464.8 | 666.0 |

2001 | Male | 10.0 | 10.0 | 964.3 | 273.0 |

2007 | Female | 10.0 | 10.0 | 195.3 | 101.5 |

2009 | Female | 10.0 | 10.0 | 224.0 | 155.7 |

2011 | Female | 10.0 | 51.0 | 679.9 | 322.5 |

2012 | Male | 10.0 | 73.5 | 343.4 | 173.5 |

2013 | Female | 10.0 | 10.0 | 229.1 | 145.1 |

2018 | Female | 10.0 | 12.3 | 243.6 | 218.8 |

2020 | Female | 10.0 | 10.0 | 73.2 | 88.0 |

2028 | Female | 10.0 | 40.8 | 814.5 | 326.0 |

2030 | Female | 10.0 | 47.0 | 286.9 | 155.1 |

2031 | Female | 10.0 | 10.0 | 166.9 | 68.8 |

2034 | Male | 10.0 | 10.0 | 474.1 | 237.0 |

2036 | Male | 10.0 | 42.9 | 306.4 | NA |

2040 | Male | 10.0 | 10.0 | 249.3 | 188.0 |

2042 | Female | 10.0 | 34.8 | 484.1 | 200.0 |

2043 | Male | 10.0 | 10.0 | 37.4 | 108.9 |

2044 | Male | 10.0 | 30.5 | 927.8 | 341.1 |

2045 | Female | 10.0 | 10.0 | 266.8 | 139.7 |

2046 | Male | 10.0 | 10.0 | 262.2 | 75.6 |

2047 | Male | 10.0 | 279.4 | 313.1 | 245.2 |

2048 | Male | 10.0 | 10.0 | 678.0 | 297.4 |

2049 | Male | 10.0 | 10.0 | 182.7 | 167.5 |

2050 | Male | 10.0 | 10.0 | 431.5 | 217.9 |

2051 | Female | 10.0 | 79.2 | 435.8 | NA |

2052 | Male | 10.0 | 10.0 | 252.1 | 151.3 |

2053 | Female | 10.0 | 27.5 | 319.3 | 184.0 |

2054 | Female | 10.0 | 35.1 | 179.2 | 117.1 |

2055 | Female | 10.0 | 28.2 | 493.1 | 342.0 |

2056 | Male | 10.0 | 33.1 | 387.3 | 466.0 |

2002 | Female | 10.0 | 77.0 | 801.1 | 617.0 |

2005 | Female | 10.0 | 19.4 | 349.3 | NA |

2006 | Female | 10.0 | 44.6 | 423.2 | 217.0 |

2008 | Female | 10.0 | 61.5 | 291.7 | 101.7 |

2010 | Female | 10.0 | 34.5 | 473.0 | NA |

2014 | Female | 10.0 | 10.0 | 15.8 | NA |

2016 | Female | 10.0 | 44.1 | 483.2 | NA |

2019 | Female | 10.0 | 10.0 | 112.0 | NA |

2022 | Female | 10.0 | 10.0 | 131.4 | 108.0 |

2023 | Female | 10.0 | 13.0 | 203.3 | NA |

2025 | Female | 58.9 | 856.3 | 2095.1 | 769.0 |

2026 | Female | 10.0 | 21.4 | 387.2 | NA |

2027 | Female | 10.0 | 10.0 | NA | 80.1 |

2003 | Male | 10.0 | 10.0 | 267.0 | NA |

2004 | Male | 10.0 | 12.6 | 1268.4 | 523.0 |

2015 | Male | 10.0 | 2876.5 | 4379.5 | 2492.0 |

2017 | Male | 10.0 | 122.9 | 587.6 | 546.5 |

2024 | Male | 10.0 | 10.0 | 10.0 | NA |

2029 | Male | 10.0 | 10.0 | 341.9 | 159.8 |

2032 | Male | 10.0 | 1143.4 | 1845.0 | 247.0 |

2033 | Male | 10.0 | 15.4 | NA | NA |

2035 | Male | 10.0 | 35.3 | 422.1 | NA |

2037 | Male | 10.0 | 26.9 | NA | NA |

2038 | Male | 10.0 | 45.9 | 391.5 | NA |

2039 | Male | 10.0 | 10.0 | 143.6 | NA |

2041 | Male | 47.8 | 2131.4 | 3044.5 | NA |

2021 | Female | 10.0 | 10.0 | 10.0 | NA |

So it appears that there are 57 unique participants. Some participants were lost to follow-up and the final date range is kind of wide (90-120 days).

One of the most important steps in any analysis is to graph the data to gain a general understanding of the trends and shape of the data.

```
dat_long %>%
ggplot(aes(day, value, group = PID))+
geom_line()+
facet_wrap(~PID)+
labs(
title = "Antibody Titres by Participate"
)
```

One things that stands out to me is participant 2015. This individual has nAbs that are much higher than the other participants.

```
dat_long %>%
mutate(tst = ifelse(PID==2015,"PID: 2015","All Others")) %>%
group_by(tst, day) %>%
summarise(mu = mean(value, na.rm=TRUE)) %>%
spread(day,mu) %>%
knitr::kable(digits = 1, caption = "Mean Value of nAbs by Group")
```

tst | 0 | 21 | 42 | 90 |
---|---|---|---|---|

All Others | 11.5 | 104.0 | 499.3 | 251.1 |

PID: 2015 | 10.0 | 2876.5 | 4379.5 | 2492.0 |

If this were me, I would go back and check those data just in case someone got a little too happy when titrating (been there). Regardless, accepting this level of variability in the individual, these kind of data lend themselves well to hierarchical modeling.

We can model the half life of the nAbs in the participants using the multilevel model framework with lme4. I will subset the data to include only the levels after the peak antibody level and include a random effect for the participant (since we know that there is within-person variability and serial dependence).

```
library(lme4)
fit_lme <- lmer(log(value) ~ day + (1|PID), data = subset(dat_long, day >21))
(half_life <- -log(2)/confint(fit_lme)[4,])
```

```
2.5 % 97.5 %
46.31701 84.74999
```

So based on this model, the half-life for nAbs is between 46.3 and 84.7 days.

```
nab_out <- function(t,hl){
2^(-t/hl)
}
sim_dat <- tibble(day = 1:365) %>%
mutate(lower_bound = nab_out(day, half_life[1]),
upper_bound = nab_out(day, half_life[2]))
sim_dat %>%
ggplot(aes(x = day))+
geom_ribbon(aes(ymin = lower_bound,
ymax = upper_bound),
fill = "orange", alpha = .5)+
geom_hline(yintercept = .25, col = "red", lty = "dashed")+
labs(
title = "Simulated Decay of nAbs",
y = "Percent of Maximum nAbs",
x = "Day After Peak nAbs"
)
```

So it appears that in roughly 3-6 months after maximum antibody levels a given vaccinated person would have 1/4 of nAb titres.

```
days_until_perc <- function(target, hl){
-log(target)*hl/log(2)
}
round(days_until_perc(.25, half_life)/30,1)
```

```
2.5 % 97.5 %
3.1 5.6
```

The big question of course is first, does this matter and in what context (protection vs infection where high levels of nAbs make sense vs symptomatic disease where a trained immune system is all that is required (thinking memory B and T cells)) and what concentration is enough at each level.

Is 5% of peak is enough then we have between 6 months and a year, perhaps in time for a teaser from infection to boost nAbs.

```
round(days_until_perc(.05, half_life)/30,1)
```

```
2.5 % 97.5 %
6.7 12.2
```

This is not my forte, but it would seem that compartmental models would be a good way to approach this problem. We would have to fix some parameters due to the data available but one could imagine that a decent data generating process could look like (at least as far as I understand how mRNA vaccines work):

- Inoculum is added to the blood stream via vaccination
- A latent phase for vaccinated individual’s cells to produce the spike protein which would be a function of the volume of inoculum and the number of cells encountered. The inoculum also has some decay function.
- A phase where the encoded protein is in the bloodstream and recognized by the immune system as a foreign body which ramps up the innate and humoral immune response.
- Eventually the spike proteins would be cleared at some rate as the immune response increases

I started down this path and then realized that I would need to make some assumptions for some of the parameters (transition rates between compartments). Being that I do not have reasonable priors for these rates, I decided not to code up this mdoel.

I read another paper on waning immunity on coronaviruses that had much more longitudinal data (Townsend et al. 2021). The key model in the paper was locked away in a proprietary Mathematica Notebook that couldn’t open without buying the software. I coded the below model based on their methods, but we will find is that we do not have enough data to fit a very good model.

```
rw("decay.stan")
```

```
data {
int<lower=1> N; // Number of Samples
int<lower=1> J; //Number of participants
vector [N] nAbs; //The concentration of the antibodies
vector [N] time; //The concentration of the antibodies
}
parameters{
real omega;
real hl;
real<lower=0> sigma;
}
transformed parameters {
vector[N] theta; //vac effect
theta[1:N] = omega + (1- omega) * exp(-hl*time[1:N]);
}
model {
hl ~ normal(0,2);
omega ~ normal(0,1);
nAbs ~ normal(theta, sigma);
}
```

```
library(cmdstanr)
mod_ex <- cmdstan_model("decay.stan")
dat_stan_base <- dat_long %>%
filter(!is.na(value)) %>%
filter(day>21)
stan_dat <- list(
N = nrow(dat_stan_base),
J = length(unique(dat_stan_base$PID)),
nAbs = dat_stan_base$value,
time = dat_stan_base$day
)
fit_ex <- mod_ex$sample(stan_dat,
refresh = 0, iter_sampling = 500,
iter_warmup = 200,
adapt_delta = .98,
max_treedepth = 15)
```

```
Running MCMC with 4 parallel chains...
Chain 1 finished in 0.3 seconds.
Chain 3 finished in 0.1 seconds.
Chain 2 finished in 0.5 seconds.
Chain 4 finished in 1.1 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.5 seconds.
Total execution time: 1.4 seconds.
```

Not a great fit as we are trying to draw a complicated curve through basically two points and there are (infinitely) many ways to do this.

```
fit_ex$summary(variables = c("omega", "hl", "sigma"))
```

```
# A tibble: 3 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 omega 0.0893 0.115 0.922 1.14 -1.06 1.38 1.28 11.5
2 hl 0.347 -0.323 1.60 1.37 -1.26 3.33 2.34 5.03
3 sigma 400. 321. 402. 476. 0.216 878. 2.32 5.03
# … with 1 more variable: ess_tail <dbl>
```

A model that is somewhat of a compromise between the compartmental model which captures more of the underlying biological process and the hierarchical model which fit the within-host variability is to model the viral decay with a direct functional form. This is similar to the approach performed by Singanayagam et al. (2021) who use this functional form for estimated RNA copies from Cycle Threshold values. The issues, as before, is that we only have a few observations per person and I don’t have strong priors, but this model is closer to what we see in the graphs where the titre increases to some point and then decay afterwards. Ideally we would have some additional measures for the participants.

This can be modeled by the following equations:

\[ v(\tau) = v_{max} \frac{(a+b)}{be^{-a(\tau-\tau_{max})}+ae^{b(\tau-\tau_{max})}} \]

The a and b parameters can be fit and influence the rate of increase and decrease of nAb tire (or virus titre for that matter).

We need to comport our data to provide the ability to identify the individual maximums and times for each participant.

And our model:

```
rw("rise.stan")
```

```
functions{
real viral_kinetics (real v_max, real a, real b, real t, real t_max){
return(v_max * (a + b) / (b * exp(-a*(t - t_max)) + a * exp(b*(t-t_max))));
}
}
data {
int<lower=1> N; // Number of records
int<lower=0> id[N]; //Participant IDs
vector<lower=0> [N] nAbs_max; // Maximum Ab concentration
vector<lower=0> [N] nAbs; //The concentration of the antibodies
vector [N] time; //The concentration of the antibodies
vector [N] t_max; //The concentration of the antibodies
real prior_a_mu;
real prior_b_mu;
real prior_a_sd;
real prior_b_sd;
int pred_window;
}
parameters {
real<lower=0> a;
real<lower=0> b;
real<lower=0> sigma;
}
transformed parameters{
vector<lower=0>[N] mu;
for( i in 1:N){
mu[i] = viral_kinetics(nAbs_max[id[i]], a, b, time[id[i]], t_max[id[i]]);
}
}
model {
a ~ normal(prior_a_mu, prior_a_sd);
b ~ normal(prior_b_mu, prior_b_sd);
sigma ~ normal(0,1);
nAbs ~ normal(mu, sigma);
}
generated quantities{
vector[pred_window] yhat;
vector[pred_window] mu_pred;
real avg_nabs_max = mean(nAbs_max);
real avg_tmax = mean(t_max);
for(i in 1:pred_window){
mu_pred[i] = viral_kinetics(avg_nabs_max, a, b, i, avg_tmax);
yhat[i] = normal_rng(mu_pred[i], sigma);
}
}
```

We can fit the model in the usual ways:

```
mod_titre <- cmdstan_model("rise.stan")
titre_dat <- list(
N = nrow(long_model_dat),
id = long_model_dat$id,
nAbs_max = long_model_dat$v_max,
nAbs = long_model_dat$value,
time= long_model_dat$day,
t_max = long_model_dat$t_max,
prior_a_mu = .8,
prior_b_mu = .1,
prior_a_sd = .05,
prior_b_sd = .05,
pred_window = 180
)
set.seed(336)
fit_titre <- mod_titre$sample(data = titre_dat,
refresh = 0, adapt_delta = .99,
max_treedepth = 15)
```

```
Running MCMC with 4 parallel chains...
Chain 1 finished in 2.2 seconds.
Chain 2 finished in 2.3 seconds.
Chain 3 finished in 2.3 seconds.
Chain 4 finished in 2.2 seconds.
All 4 chains finished successfully.
Mean chain execution time: 2.3 seconds.
Total execution time: 2.4 seconds.
```

```
fit_titre$summary(variables = c("a", "b", "sigma"))
```

```
# A tibble: 3 × 10
variable mean median sd mad q5 q95 rhat ess_bulk
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 a 0.801 0.801 0.0511 0.0501 0.714 0.885 1.00 2165.
2 b 0.102 0.0996 0.0479 0.0505 0.0252 0.182 1.00 1447.
3 sigma 91.9 91.9 0.499 0.494 91.1 92.7 1.01 2466.
# … with 1 more variable: ess_tail <dbl>
```

Note that there is a fair bit of variance in our model represented by zero. Additionally, we have some values of the predicted titre below zero. We could resolve this by changing the distribution to use a distribution that is always positive like a lognormal distribution. This would require some careful transformation of our data and associated parameters. We’ll stick with this formulation for now.

We can now visualize our results from this model.

```
library(posterior)
library(ggdist)
pred_nab <- tidybayes::gather_draws(fit_titre$draws(variables = "yhat"),yhat[i] ) %>%
dplyr::ungroup() %>%
dplyr::group_by(i) %>%
ggdist::median_qi(.value, .width = c(.5, .8, .95))
pred_nab %>%
ggplot(aes(i, .value))+
geom_lineribbon(aes(ymin = .lower, ymax = .upper))+
scale_fill_brewer()+
labs(
title = "Estimated Quantity of nAbs with Viral Evolution Model",
y = "Quantity",
x = "Days",
fill = "Proportion of Scenarios"
)+
geom_hline(yintercept = 0, col = "orange", lty = "dashed")+
theme(legend.position = "top")
```

Again, we see evidence of detectable tire above 180 days, but there is a decent bit of variability.

If we had additional measurements per person and some additional information on the kinetics, In the time it took me to write this post several other papers have come out with much more information that would be fun to analyze that find that immunity is likely long lasting (Mateus et al., n.d.; Grifoni et al. 2021).

Corbett, Kizzmekia S., Martha C. Nason, Britta Flach, Matthew Gagne, Sarah O’Connell, Timothy S. Johnston, Shruti N. Shah, et al. 2021. “Immune Correlates of Protection by mRNA-1273 Vaccine Against SARS-CoV-2 in Nonhuman Primates.” *Science*, July. https://doi.org/10.1126/science.abj0299.

Grifoni, Alba, John Sidney, Randi Vita, Bjoern Peters, Shane Crotty, Daniela Weiskopf, and Alessandro Sette. 2021. “SARS-CoV-2 Human t Cell Epitopes: Adaptive Immune Response Against COVID-19.” *Cell Host & Microbe* 29 (7): 1076–92. https://doi.org/10.1016/j.chom.2021.05.010.

Mateus, Jose, Jennifer M. Dan, Zeli Zhang, Carolyn Rydyznski Moderbacher, Marshall Lammers, Benjamin Goodwin, Alessandro Sette, Shane Crotty, and Daniela Weiskopf. n.d. “Low-Dose mRNA-1273 COVID-19 Vaccine Generates Durable Memory Enhanced by Cross-Reactive t Cells.” *Science* 0 (0): eabj9853. https://doi.org/10.1126/science.abj9853.

Singanayagam, Anika, Seran Hakki, Jake Dunning, Kieran J. Madon, Michael A. Crone, Aleksandra Koycheva, Nieves Derqui-Fernandez, et al. 2021. “Community Transmission and Viral Load Kinetics of the SARS-CoV-2 Delta (b.1.617.2) Variant in Vaccinated and Unvaccinated Individuals in the UK: A Prospective, Longitudinal, Cohort Study.” *The Lancet Infectious Diseases* 0 (0). https://doi.org/10.1016/S1473-3099(21)00648-4.

Suthar, Mehul S., Prabhu S. Arunachalam, Mengyun Hu, Noah Reis, Meera Trisal, Olivia Raeber, Sharon Chinthrajah, et al. 2021. “Durability of Immune Responses to the BNT162b2 mRNA Vaccine.” https://doi.org/10.1101/2021.09.30.462488.

Townsend, Jeffrey P., Hayley B. Hassler, Zheng Wang, Sayaka Miura, Jaiveer Singh, Sudhir Kumar, Nancy H. Ruddle, Alison P. Galvani, and Alex Dornburg. 2021. “The Durability of Immunity Against Reinfection by SARS-CoV-2: A Comparative Evolutionary Study.” *The Lancet Microbe* 0 (0). https://doi.org/10.1016/S2666-5247(21)00219-6.

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Text and figures are licensed under Creative Commons Attribution CC BY 4.0. Source code is available at https://github.com/medewitt/medewitt.github.io, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

For attribution, please cite this work as

DeWitt (2021, Nov. 11). Michael DeWitt: Neutralizing Antibody Titres. Retrieved from https://michaeldewittjr.com/programming/2021-10-11-neutralizing-antibody-titres/

BibTeX citation

@misc{dewitt2021neutralizing, author = {DeWitt, Michael}, title = {Michael DeWitt: Neutralizing Antibody Titres}, url = {https://michaeldewittjr.com/programming/2021-10-11-neutralizing-antibody-titres/}, year = {2021} }