Neutralizing Antibody Titres

R antibodies nAbs Covid-19 Immunology Bayes Stan

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.

Michael DeWitt https://michaeldewittjr.com
2021-11-11

Reading the Pre-Prints

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)

So Look at the Data

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"
  )
Antibody Titres Over Time by Participant

Figure 1: Antibody Titres Over Time by Participant

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")
Table 1: 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.

So What’s the Half Life?

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 

What about compartmental models?

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):

  1. Inoculum is added to the blood stream via vaccination
  2. 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.
  3. 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.
  4. 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.

Other Method

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>

Viral Titre?

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

nab_kinetics <- function(v_max, a,b,t,tmax){
  v_max * (a+b)/(b*exp(-a*(t-tmax))+a*exp(b*(t-tmax)))
}

plot(nab_kinetics(50, a = .8, b= .1, t = 1:70, tmax = 20), type = "l",
     ylab = "nAb Titre", main = "Simulated Results")

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

long_model_dat <- dat_long %>% 
  filter(!is.na(value)) %>% 
  group_by(PID) %>% 
  mutate(v_max = max(value)) %>% 
  mutate(t_max =  max((v_max == value)*day))%>% 
  ungroup() %>% 
  mutate(id = as.integer(as.factor(PID)))

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.

Conclusions

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.

References

Corrections

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

Reuse

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

Citation

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}
}