In mrc-ide/safir: Individual-Based Model of COVID Transmission

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

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

Antibody dynamics

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

Titre associated with dose or infection

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.

```{R, eval=FALSE} 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.

```{R, eval=FALSE}
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

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:

```{R, eval=FALSE}

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:

```{R, eval=FALSE}
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

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.

```{R, eval=FALSE}

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:

```{R, eval=FALSE}
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.

Efficacy against onward infection

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

```{R, eval = FALSE}

group infectious persons by age

inf_ages <- get_inf_ages(infection_bset = infectious, variables = variables, parameters = parameters)

The relevant (R) code is here:

```{R, eval=FALSE}
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.

mrc-ide/safir documentation built on Aug. 2, 2022, 10:47 a.m.