Simulating a clinical study with a mechanistic model in the background

Description of virtual study

• Increased doses for leukemia in comparison to “standard” doses with some substance X is to be studied
• Randomized, controlled, double-blind clinical trial
• Groups:
  • standard regimen: X 1 mg - 0 - 2 mg - 0 for 10 days
  • adjusted regimen: X 2 mg - 0 - 3 mg - 0 for 10 days
• Leukemic burden is assumed to be lognormally distributed initially
• Postleukemic burden is assumed to be: Leukemic burden pretherapy divided by area under the plasma level curve
• Some plausible pharmacokinetic model is assumed, with rate constants depending on age and weight
• Primary endpoint: Molecular remission (defined as postleukemic burden < 0.0012)

Definition of study population

set.seed(123)
nr_of_patients <- 200
patients <- tibble(id = 1:nr_of_patients,
                   group = rep(c('standard_dosing', 'adjusted_dosing'), each = nr_of_patients / 2),
                   age = rnorm(nr_of_patients, mean = 60, sd = 10),
                   weight = rnorm(nr_of_patients, mean = 80, sd = 10),
                   age_scaled = scale(age)[, 1],
                   weight_scaled = scale(weight)[, 1],
                   sex = sample(c('male', 'female'), 
                                size = nr_of_patients, 
                                replace = TRUE),
                   leukemic_burden_pretherapy = rlnorm(nr_of_patients))

patients |> 
  sample_n(5) |> 
  kable(align = 'c')

id	group	age	weight	age_scaled	weight_scaled	sex	leukemic_burden_pretherapy
182	adjusted_dosing	72.63185	81.90230	1.3483988	0.1487158	female	2.670092
13	standard_dosing	64.00771	92.32476	0.4340111	1.1952301	male	2.416850
68	standard_dosing	60.53004	62.43473	0.0652855	-1.8060152	female	1.798109
127	adjusted_dosing	62.35387	79.39178	0.2586592	-0.1033647	female	1.186002
25	standard_dosing	53.74961	76.36343	-0.6536207	-0.4074403	male	3.048497

Definition of pharmacokinetic model

calc_plasma_levels <- function(applications, patient, t_max = 10) {
  plasma <- tibble(timepoint = 1:t_max, 
                   level = 0, 
                   id = patient$id, 
                   group = patient$group) 
  for(i in 1:nrow(applications)) {
    ii <- applications$time[i]
    plasma$level[ii:t_max] <- plasma$level[ii:t_max] + applications$dose[i]
  }
  degraded <- 0
  for(i in 2:t_max) {
    degraded <- degraded + 
      plasma$level[i - 1] * 
      (0.05 + (patient$age_scaled + patient$weight_scaled) / 100)
    plasma$level[i] <- plasma$level[i] - degraded
  }
  plasma
}

Definition of dosing schemes

standard_dosing <- tibble(time = c(12 * 0:19 + 1), 
                            dose = rep(c(1, 2), length.out = 20))
adjusted_dosing <- tibble(time = c(12 * 0:19 + 1), 
                            dose = rep(c(2, 3), length.out = 20))

Calculation of plasma levels

plasma_levels <- tibble()
for(i in 1:nr_of_patients){
  if(patients$group[i] == 'standard_dosing') {
    current_dosing <- standard_dosing
  } else {
    current_dosing <- adjusted_dosing
  }
  plasma_levels <- plasma_levels |> 
    bind_rows(calc_plasma_levels(
      applications = current_dosing, 
      patient = patients[i, ],
      t_max = 300))
}

Analysis

plasma_levels |> 
  filter(id %in% sample(nr_of_patients, 20)) |> 
  ggplot(aes(timepoint, level, group = factor(id), colour = group)) + geom_line() + 
  labs(title = 'Plasma levels of randomly chosen individuals', colour = 'Patient id')

patients <- patients |> 
  full_join(group_by(plasma_levels, id) |> summarize(auc = sum(level)))

Joining with `by = join_by(id)`

patients <- patients |> 
  mutate(leukemic_burden_posttherapy = leukemic_burden_pretherapy / auc,
         molecular_remission = leukemic_burden_posttherapy > 1.2e-3)

summary(patients$leukemic_burden_posttherapy)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
8.897e-05 5.727e-04 1.214e-03 1.981e-03 2.403e-03 2.008e-02 

# Contigency table of group vs. molecular_remission (yes / no)
(study_outcome <- xtabs(~ group + molecular_remission, data = patients))

                 molecular_remission
group             FALSE TRUE
  adjusted_dosing    61   39
  standard_dosing    38   62

# Calculate study p-value
chisq.test(study_outcome)

    Pearson's Chi-squared test with Yates' continuity correction

data:  study_outcome
X-squared = 9.681, df = 1, p-value = 0.001862

Conclusion

• Here, a hypothetical clinical study has been simulated
• a lot of fine-tuning and a lot of knowledge about the underlying mechanisms is needed
• hard to generalize that for the next application - probably, a similar amount of fine-tuning will be needed