Busulfan MIPD
busulfan_mipd.Rmd
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)