Busulfan MIPD • mipdtrial

In this example script, we will estimate target attainment for busulfan.

To account for model misspecification, we will use one model to simulate (McCune, 2014) and one model to estimate (Shukla, 2020). This is a simplified version of a published study in which we compared non-compartmental analysis (NCA) with MAP Bayesian estimation for AUC-guided busulfan dosing.

library(mipdtrial)
if(!require(pkbusulfanmccune)) {
  PKPDsim::install_default_literature_model("pk_busulfan_mccune")
  library(pkbusulfanmccune)
}
if (!require(pkbusulfanshukla)) {
  PKPDsim::install_default_literature_model("pk_busulfan_shukla")
  library(pkbusulfanshukla)
}

# For example data
library(NHANES)

# For data handling/plotting
library(dplyr)
library(tidyr)
library(ggplot2)

1. Define trial design

In this simple trial, patients will be dosed once per day for four days, with 4 levels collected per day, at 3.5, 4, 6 and 8 hours post-dose. We will adjust doses 2, 3, and 4 to achieve a target cumulative AUC of 90 by t = 192 hours.

tdm_design <- create_sampling_design(
  offset = rep(c(3.5, 4, 6, 8), 4),
  at = rep(1:4, each = 4),
  anchor = "dose"
)
target_design <- create_target_design(
  targettype = "cum_auc", 
  targetvalue = 90 * 1000,       # unit conversion
  time = 192
)
update_design <- create_regimen_update_design(
  at = c(2, 3, 4),
  anchor = "dose",
  dose_optimization_method = map_adjust_dose
)

2. Create a set of digital patient covariates

We will use the NHANES data set for covariate distributions. This is a data set of five thousand individuals collected by the US National Center for Health Statistics. See package documentation for more details.

We will also use the dplyr package for some basic data processing. Since every data set is un-tidy in its own way (“Tidy datasets are all alike, but every messy dataset is messy in its own way.” –– Hadley Wickham), data cleaning should be handled outside of the mipdtrial tools. The functions expect all patient covariates to be numeric.

dat <- NHANES::NHANES %>%
  # take just the first record per patient
  group_by(ID) %>%
  slice(1) %>%
  ungroup() %>%
  # convert columns to numeric
  mutate(Gender = ifelse(Gender == "male", 1, 0)) %>% # as defined in models
  select(ID, Gender, Age, Weight, Height) %>%
  # include only patients with all information available
  filter(!is.na(Gender), !is.na(Weight), !is.na(Height), !is.na(Age))

The next step is to make sure we have all the covariates we need. Busulfan has a known time-dependent relationship with clearance, and so the McCune model takes a covariate that indicates at what time this relationship should start (T_CL_EFF; useful for test doses that should not impact busulfan clearance). We will set this to zero for all patients. The Shukla model also takes conditioning regimen as a covariate, since co-medication with Clofarabine and Fludarabine was associated with a change in clearance (REGI). We will also set this to zero for all patients.

Other models might require fat-free mass or other calculated covariates. This would be a good time to do that sort of processing on your data set!

You may also need to convert between units (height in m to cm, for example.)

dat$T_CL_EFF <- 0
dat$REGI <- 0

Here are the first few rows of our data set:

head(dat)
#> # A tibble: 6 × 7
#>      ID Gender   Age Weight Height T_CL_EFF  REGI
#>   <int>  <dbl> <int>  <dbl>  <dbl>    <dbl> <dbl>
#> 1 51624      1    34   87.4   165.        0     0
#> 2 51625      1     4   17     105.        0     0
#> 3 51630      0    49   86.7   168.        0     0
#> 4 51638      1     9   29.8   133.        0     0
#> 5 51646      1     8   35.2   131.        0     0
#> 6 51647      0    45   75.7   167.        0     0

We will also only look at the first 10 patients in our data set, in the interest of speed.

ids <- dat$ID[1:10]

3. Load model definitions

We need to define the model parameters and variability terms for our simulation. For this study, we will use the ones defined in the models.

# Set up simulation model ("true")
model_sim <- pkbusulfanmccune::model()
parameters_sim <- pkbusulfanmccune::parameters()
omega_sim <- pkbusulfanmccune::omega_matrix()
ruv_sim <- pkbusulfanmccune::ruv()

# Set up estimation model (model for MIPD)
model_est <- pkbusulfanshukla::model()
parameters_est <- pkbusulfanshukla::parameters()
omega_est <- pkbusulfanshukla::omega_matrix()
ruv_est <- pkbusulfanshukla::ruv()

We also need to link the covariates in our data set to the covariates expected in the model.

# check which covariates are required for your model using 
# `attr(model_sim, "covariates")`

# check which covariates are in your data set using `colnames(dat)`
cov_mapping <- c(
  AGE = "Age", 
  WT = "Weight", 
  HT = "Height", 
  SEX = "Gender",
  T_CL_EFF = "T_CL_EFF",
  REGI = "REGI"
)

4. Initialize output metrics of interest

We want to record both “true” cumulative AUC and estimated cumulative AUC, along with patient identifiers (i.e., id). We will pre-populate a data frame rather than iteratively binding rows, since that is more performant.

results <- data.frame(
  id = ids,
  iter = 1,
  true_auc = NA_real_,
  est_auc = NA_real_
)

5. Simulate a trial!

Each patient will start at a dose of 3.2 mg/kg body weight. We will then use mipdtrial::dose_grid_search to optimize doses to attain a target of 90 mg*h/L.

For inter-individual variability terms, we will use the inter-individual variability described by our “true” model, since this should reflect the “true” distribution in patient pharmacokinetics. Later, we will estimate the individual PK parameters using our estimation model and the measured TDMs we “collect”.

For residual error terms, we will use the residual error of our “true” model, since this is supposed to reflect the “true” unexplained error (assay error, etc.). This function creates a data frame of error to add to each “true” simulated concentration to produced a measured therapeutic drug monitoring sample.

set.seed(1) # important for reproducibility
for (i in ids) {
  # get patient covariates
  covs <- create_cov_object(
    dat[dat$ID == i,],
    mapping = cov_mapping
  )
  # create initial dosing regimen (see PKPDsim::new_regimen for more info)
  reg <- PKPDsim::new_regimen(
    amt = 3.2 * covs$WT$value, # we will update this each "day" of the trial.
    interval = 24,
    n = 4,
    t_inf = 3,
    type = "infusion"
  )
  # randomly draw individual PK parameters
  pars_true_i <- generate_iiv(
    sim_model = model_sim,
    omega = omega_sim,
    parameters = parameters_sim
  )
  
  treatment_summary <- sample_and_adjust_by_dose(
    sampling_design = tdm_design,
    regimen_update_design = update_design,
    target_design = target_design,
    regimen = reg,
    covariates = covs,
    pars_true_i = pars_true_i,
    sim_model = model_sim,
    sim_ruv = ruv_sim,
    est_model = model_est,
    parameters = parameters_est,
    omega = omega_est,
    ruv = ruv_est
  )
  # extract metrics of interest
  results$true_auc[results$id == i] <- calc_auc_from_regimen(
    regimen = treatment_summary$final_regimen,
    parameters = pars_true_i, # true patient parameters
    model = model_sim,
    target_design = target_design,
    covariates = covs
  )
  final_parameters <- treatment_summary$additional_info[[length(treatment_summary$additional_info)]]
  results$est_auc[results$id == i] <- calc_auc_from_regimen(
    regimen = treatment_summary$final_regimen,
    parameters = final_parameters, # estimated patient PK
    model = model_est,
    target_design = target_design,
    covariates = covs
  )
}

Here are the first few rows of our simulation results:

head(results)
#>      id iter  true_auc  est_auc
#> 1 51624    1  83063.85 90006.79
#> 2 51625    1 103788.19 89998.38
#> 3 51630    1  83271.64 90045.76
#> 4 51638    1  95389.50 89999.14
#> 5 51646    1  82884.90 89965.09
#> 6 51647    1 104570.85 89964.21

6. Analyze results

How well did our patients get to target? It looks like we estimated that target attainment would be very high (all patients have an AUC very close to 90 mg*h/L), but “true” AUC was higher than that due to some model misspecification!

results %>%
  pivot_longer(
    c(true_auc, est_auc),
    names_to = "auc_type",
    values_to = "auc"
  ) %>%
  ggplot() +
    aes(x = auc_type, y = auc/1000) +
    geom_boxplot() +
    geom_hline(yintercept = 90, color = "navyblue", linetype = "dashed") +
    scale_y_continuous(limit = c(0, NA)) +
    theme_minimal() +
    labs(
      y = "Cumulative AUC (mg\u00B7h/L)",
      x = "AUC type"
    )

Even with some model misspecification, target attainment (within 20% of target AUC) was still high:

target_attainment <- results %>%
  mutate(
    ontarget = is_on_target(true_auc, target_design)) %>%
  summarize(proportion_ontarget = mean(ontarget)) %>%
  pull(proportion_ontarget)
target_attainment <- paste0(round(100 * target_attainment), "%")

Overall, 90% of patients had a “true” AUC of 90 mg*h/L.

For an in-depth analysis of how differences between these two models impact target attainment, see Hughes et al., J PKPD (2024)