vignettes/antibody_efficacy.Rmd
antibody_efficacy.RmdThis vignette will explain how antibody titre dynamics are modeled in safir, and how they are affected by vaccination and (optionally) natural immunity from infection bouts. Furthermore it will explain how protective efficacy against infection and severe disease, and the effect on infectiousness is calculated from an individual’s antibody titre.
The variables needed to store antibody titre are created in create_vaccine_variables, in R/variables_vaccines.R, attached to the variable list variables$ab_titre. The antibody titre decays and that process is vaccine_ab_titre_process (R/efficacy_vaccination.R) or natural_immunity_ab_titre_process (R/ efficacy_naturalimmunity.R), depending on if you are modeling natural immunity or not.
Each day, an individual’s antibody titre decays at a rate given by the vector parameters$dr_vec, giving the daily rate of decay since last vaccine dose (or recovery from infection bout, for natural immunity). If the person’s time since last dose (or infection) is greater than the length of the vector, the last value is used. The vector is set in get_vaccine_ab_titre_parameters (R/parameters_vaccine.R).
The only way that antibody titre can increase is through a vaccine dose or, optionally, natural infection. For vaccination, the mean titre associated with each dose is associated with the particular vaccine product, and is chosen in get_vaccine_ab_titre_parameters. At this point we note that there are two ways to model the actual titre drawn for each individual receiving a dose. By default the value in mu_ab_list for the chosen vaccine product is used. However if correlated = TRUE is used, then the individual’s antibody titre sampled from the previous dose is multiplied by the ratio of the current to previous dose mean titre values.
The function which samples antibody titre for each dose is schedule_dose_vaccine and is in R/events_vaccination.R. Note that if correlated = TRUE was selected then there is an additional variable zdose which holds the previous antibody titre associated with the last vaccine dose for each individual.
Using correlated doses increases memory requirements for the model, and produces nearly identical population level trajectories of antibody titre as uncorrelated, so its use should be weighed against other computational factors.
schedule_dose_vaccine <- function(timestep, variables, target, dose, parameters) {
variables$dose_num$queue_update(value = dose,index = target)
variables$dose_time$queue_update(values = timestep, index = target)
if (inherits(target,"Bitset")) {
n <- target$size()
} else {
n <- length(target)
}
if (parameters$correlated) {
if (dose > 1) {
# correlated doses > 1; use ratio of mean titre
zdose_prev <- variables$zdose$get_values(index = target)
zdose_prev <- log10(exp(zdose_prev)) # transform back to log scale
zdose <- log(10^zdose_prev * (parameters$mu_ab[dose] / parameters$mu_ab[dose-1]))
variables$zdose$queue_update(values = zdose, index = target)
} else {
# initial dose
zdose <- log(10^rnorm(n = n, mean = log10(parameters$mu_ab[dose]),sd = parameters$std10))
variables$zdose$queue_update(values = zdose, index = target)
}
} else {
# uncorrelated doses (also use for dose 1 of correlated dose titre)
zdose <- log(10^rnorm(n = n, mean = log10(parameters$mu_ab[dose]),sd = parameters$std10))
}
variables$ab_titre$queue_update(values = zdose, index = target)
}For antibody titre driven by natural immunity, the recovered (R) class of completely immune individuals used in the compartmental model to simulated complete immunity with Exponential or Erlang distributed wait time before transition back to susceptible (S) is now simply used as a state from which to compute the antibody titre associated with that individual’s natural immune response to infection.
This requires two additional variable objects, num_inf which tracks how many infections each individual has experienced, and inf_time, the time since last infection. These variables are made in create_natural_immunity_variables (R/variables_naturalimmunity.R).
Upon entering the R compartment, individuals will transition back to S on the next time step (as time steps are usually much less than one day, this is not a problem). During that time event listeners are called which draw the antibody titre associated with this infection bout and update ab_titre. Of note is that the vector parameters$mu_ab_infection controls the mean antibody titre associated with each re-infection event, in case antibody response varies from the first to subsequent re-infections, and can be set by the user. The user can optionally specify parameters$std10_infection to control the variance of the sampled values, but if not specified it will use the std10 parameter for vaccines.
Antibody titre boost associated with recovery can either be overwriting or additive, controlled by the additive parameter in attach_event_listeners_natural_immunity.
The function attach_event_listeners_natural_immunity() is found in R/events_naturalimmunity.R.
attach_event_listeners_natural_immunity <- function(variables, events, parameters, dt, additive = FALSE) {
stopifnot(is.logical(additive))
stopifnot(!is.null(parameters$mu_ab_infection))
stopifnot(is.finite(parameters$mu_ab_infection))
if (is.null(parameters$std10_infection)) {
std10_infection <- parameters$std10
} else {
stopifnot(length(parameters$std10_infection) == 1L)
stopifnot(is.finite(parameters$std10_infection))
std10_infection <- parameters$std10_infection
}
# recovery: handle 1 timestep R->S and update ab titre for immune response
if (length(events$recovery$.listeners) == 2) {
events$recovery$.listeners[[2]] <- NULL
}
# they go from R to S in 1 time step
events$recovery$add_listener(
function(timestep, target) {
events$immunity_loss$schedule(target = target, delay = rep(1, target$size()))
}
)
# effect of infection on NAT
if (additive) {
events$recovery$add_listener(
function(timestep, target) {
# update inf_num
inf <- variables$inf_num$get_values(target) + 1L
variables$inf_num$queue_update(values = inf, index = target)
# update last time of infection
variables$inf_time$queue_update(values = timestep, index = target)
# get NAT values and convert to linear scale
current_ab_titre <- variables$ab_titre$get_values(index = target)
current_ab_titre <- exp(current_ab_titre)
# draw NAT boost on linear scale
inf[inf > length(parameters$mu_ab_infection)] <- length(parameters$mu_ab_infection)
zdose <- 10^rnorm(n = target$size(), mean = log10(parameters$mu_ab_infection[inf]),sd = std10_infection)
new_ab_titre <- current_ab_titre + zdose
# back to ln scale, and impose max value constraint
new_ab_titre <- log(new_ab_titre)
new_ab_titre <- pmin(new_ab_titre, parameters$max_ab)
# queue NAT update
variables$ab_titre$queue_update(values = new_ab_titre, index = target)
}
)
} else {
events$recovery$add_listener(
function(timestep, target) {
# update inf_num
inf <- variables$inf_num$get_values(target) + 1L
variables$inf_num$queue_update(values = inf, index = target)
# draw ab titre value
inf[inf > length(parameters$mu_ab_infection)] <- length(parameters$mu_ab_infection)
zdose <- log(10^rnorm(n = target$size(), mean = log10(parameters$mu_ab_infection[inf]),sd = std10_infection))
zdose <- pmin(zdose, parameters$max_ab)
variables$ab_titre$queue_update(values = zdose, index = target)
# update last time of infection
variables$inf_time$queue_update(values = timestep, index = target)
}
)
}
}Efficacy against infection is used to adjust the per-capita force of infection in the function infection_process_vaccine (R/process_infection_vaccine.R). The relevant code snippet is here:
# get infection modifier and ages
ab_titre <- variables$ab_titre$get_values(susceptible)
infection_efficacy <- vaccine_efficacy_infection_cpp(ab_titre = ab_titre,parameters = parameters)
ages <- variables$discrete_age$get_values(susceptible)
# FoI for each susceptible based on their age group
lambda <- lambda + (lambda_age[ages] * infection_efficacy)The function that computes efficacy against infection is vaccine_efficacy_infection (R/efficacy_vaccination.R) but for speed we use vaccine_efficacy_infection_cpp (src/efficacy_vaccination.cpp).
The relevant (R) code is here:
nt <- exp(ab_titre[finite_ab])
ef_infection[finite_ab] <- 1 / (1 + exp(-parameters$k * (log10(nt) - log10(parameters$ab_50)))) # reported efficacy in trials
ef_infection[finite_ab] <- 1 - ef_infection[finite_ab]nt is the antibody titre of each person converted back to the linear scale. The trial efficacy is \(\frac{1}{1 + e^{-k(\log_{10}(nt) - \log_{10}(ab_{50}))}}\). The protective efficacy is \(1\) minus this.
Note that we only compute this for individuals with finite antibody titre, safir uses -Inf as a null placeholder for persons who do not yet have antibodies.
Efficacy against severe disease adjusts the probability of a person once infected to experience a severe disease outcome. It is used in create_exposure_scheduler_listener_vaccine (R/events_exposure_vaccine.R), and the relevant code snippet is here. Note that we first need to compute efficacy against infection, because efficacy against severe disease is conditional upon this.
# vaccine efficacy against severe disease
ab_titre <- variables$ab_titre$get_values(hosp)
infection_efficacy <- vaccine_efficacy_infection_cpp(ab_titre = ab_titre,parameters = parameters)
severe_efficacy <- vaccine_efficacy_severe_cpp(ab_titre = ab_titre,ef_infection = infection_efficacy,parameters = parameters)
# sample those with severe disease
hosp$sample(prob_hosp * severe_efficacy)The function that computes efficacy against infection is vaccine_efficacy_severe (R/efficacy_vaccination.R) but for speed we use vaccine_efficacy_severe_cpp (src/efficacy_vaccination.cpp).
The relevant (R) code is here:
finite_ab <- which(is.finite(ab_titre))
nt <- exp(ab_titre[finite_ab])
ef_severe_uncond <- 1 / (1 + exp(-parameters$k * (log10(nt) - log10(parameters$ab_50_severe))))
ef_severe[finite_ab] <- 1 - ((1 - ef_severe_uncond)/ef_infection[finite_ab]) # 1 - (1 - ef_infection) goes from hazard reduction scale to efficacy, simplifies to ef_infection
ef_severe[finite_ab] <- 1 - ef_severe[finite_ab]nt is the linear scale antibody titre. Unconditional protection against severe disease is computed as \(\frac{1}{1 + e^{-k(\log_{10}(nt) - \log_{10}(ab\:\text{severe}_{50}))}}\). Because this is computing for individuals who have been infected we need to make it conditional on being infected, which we do by 1 - ((1 - ef_severe_uncond)/(1 - (1 - ef_infection))) where 1 - ef_infection goes from hazard reduction scale to trial efficiency scale, and it simplifies to ef_infection.
Finally we set ef_severe <- 1 - ef_severe so we are on the proper hazard reduction scale for efficacy against severe disease.
Antibody titre can (optionally) be used to adjust the contribution of each infectious person to the overall force of infection on susceptible persons. It can be turned on or off, and strength of the effect controlled with parameters set in get_vaccine_ab_titre_parameters.
It affects the computation of the number of infectious persons by age group in infection_process_vaccine (R/process_infection_vaccine.R):
# group infectious persons by age
inf_ages <- get_inf_ages(infection_bset = infectious, variables = variables, parameters = parameters)The relevant (R) code is here:
nt <- exp(ab_titre[finite_ab])
ef_transmission[finite_ab] <- 1 / (1 + exp(-parameters$k * (log10(nt / parameters$nt_transmission_factor) - log10(parameters$ab_50)))) # reported efficacy in trials
ef_transmission[finite_ab] <- 1 - ef_transmission[finite_ab]nt is the antibody titre of each person converted back to the linear scale. The trial efficacy is \(\frac{1}{1 + e^{-k(\log_{10}(nt/\text{nt_transmission_factor}) - \log_{10}(ab_{50}))}}\). The effect on each person’s onward contribution of infectiousness is 1 - ef_transmission.