Add compartments cmt to regimens in nm_to_regimen

[dominic-irx]

This PR changes nm_to_regimen so that it uses a map for compartment to dose type. This is useful for mipdeval, where this is needed.

This breaks backward compatibility because we now want an explicit declaration of the type of dosing instead of inferring from presence/absence of columns.

This function is used in two places:

1. PKPDposterior::nonmem_to_stan_data:40, which we haven't touched in a while
2. mipdeval::parse_input_data:45, which is where we have the temporary fix

[jasmineirx]

It will be great to have this functionality available, and the testing looks good for this approach. I have one question about the design choice.

I see you highlighted this as a breaking change. Is there a way it could be elegantly designed such that it is not a breaking change? For example, in new_regimen, if type and t_inf are not specified, we assume the doses are administered with "type" = "bolus" and "t_inf" = 0. If cmt_mapping is not provided and cmt is not specified, in apply_duration_scale we assign the dose to the first compartment. Previous behavior for this function was to assume that the type was bolus and that if rate/t_inf was non-zero then type was infusion. It seems to me like we could (without too much undue complexity) implement a form of this logic where, if CMT is not present in the data, we use RATE to determine if the dose is bolus or infusion (with cmt being 1 in that case). This would be consistent with logic elsewhere in the package and within this function.

It might also be a good idea to move cmt_mapping to be the final argument so that old pass-by-position argument calls are not broken.

[jasmineirx]

in general lgtm, a couple suggestions for clean-up. what about a case where dose_cmts is defined but the data doesn't have CMT?

[jasmineirx]

you might also want to add a check for !is.null(tmp$cmt) since the code below uses it.

[jasmineirx]

I think we could skip this argument. if it is null, PKPDsim::sim will use the compartment definitions of the model file, which is probably safer.

[dominic-irx]

I can see this leading to a situation where you specify CMT = 1 but PKPDsim keeps giving you CMT = 2 and you have no idea where it is coming from. Explicit is better than implicit, and this function is basically not used anywhere except in PKPDposterior (which is essentially unmaintained) and in mipdeval, where the original fix would have already worked.

[jasmineirx]

Could you help me understand this case with an example? It is possible that NONMEM-style data and PKPDsim model definitions use different compartments for one dose type.

[dominic-irx]

https://github.com/InsightRX/mipdeval/pull/25/files

This is the only relevant place this function is currently being used, where we want to specify which compartments are from which doses using the NONMEM data set. Leaving it to be assumed by the model downstream, especially when CMT should be specified in the data set anyway, is not optimal imo.

[dominic-irx]

I have completely overhauled this PR. As discussed with @roninsightrx , we are going to rely on CMT being passed into PKPDsim::regimen() and having that be interpreted correctly by the model.

First, I had to check that passing compartment and type in new_regimen will result in compartment being respected by the simulator. For voriconazole Friberg, oral is compartment 1 and IV is compartment 2. The plot shows the cmt is ultimately what matters for the simulation.

library(PKPDsim)
library(pkvoriconazolefriberg)
library(ggplot2)
library(dplyr)
reg_iv <- new_regimen(
  amt = 400,
  interval = 12,
  cmt = 2,
  t_inf = 1
)

reg_oral <- new_regimen(
  amt = 400,
  interval = 12,
  cmt = 1
)

mod <- pkvoriconazolefriberg::model()
par <- pkvoriconazolefriberg::parameters()
omega <- pkvoriconazolefriberg::omega_matrix()

dat_iv <- sim(
  seed = 604,
  ode = mod,
  parameters = par,
  omega = omega,
  covariates = list(
    "WT" = 70,
    "AGE" = 14,
    "CYP2C19unknown" = 0,
    "CYP2C19a1a2" = 0,
    "CYP2C19a2a2" = 0,
    "CYP2C19a2a3" = 0,
    "CYP2C19a1a3" = 0,
    "CYP2C19a3a3" = 0
  ),
  n_ind = 1,
  regimen = reg_iv,
  t_obs = seq(0, 36, 0.1)
)

dat_po <- sim(
  seed = 604,
  ode = mod,
  parameters = par,
  omega = omega,
  covariates = list(
    "WT" = 70,
    "AGE" = 14,
    "CYP2C19unknown" = 0,
    "CYP2C19a1a2" = 0,
    "CYP2C19a2a2" = 0,
    "CYP2C19a2a3" = 0,
    "CYP2C19a1a3" = 0,
    "CYP2C19a3a3" = 0
  ),
  n_ind = 1,
  regimen = reg_oral,
  t_obs = seq(0, 36, 0.1)
)

tmp <- data.frame(
  amt = c(400, 400, 400),
  times = c(0, 12, 24),
  cmt = c(1, 1, 2),
  rate = c(0, 0, 400)
)
reg_both <- new_regimen(
  amt = tmp$amt,
  times = tmp$time,
  cmt = tmp$cmt,
  t_inf = ifelse(tmp$rate == 0, 0, tmp$amt/tmp$rate)
)

dat_both <- sim(
  seed = 604,
  ode = mod,
  parameters = par,
  omega = omega,
  covariates = list(
    "WT" = 70,
    "AGE" = 14,
    "CYP2C19unknown" = 0,
    "CYP2C19a1a2" = 0,
    "CYP2C19a2a2" = 0,
    "CYP2C19a2a3" = 0,
    "CYP2C19a1a3" = 0,
    "CYP2C19a3a3" = 0
  ),
  n_ind = 1,
  regimen = reg_both,
  t_obs = seq(0, 36, 0.1)
)

dat_iv |>
  mutate(type = "IV") |>
  bind_rows(dat_po |> mutate(type = "PO")) |>
  bind_rows(dat_both |> mutate(type = "Both")) |>
  filter(comp == "obs") |>
  ggplot(aes(x = t, y = y, color = type)) +
  geom_point() +
  geom_line()

Next, I checked whether this function nm_to_regimen works in the same way.

# NM to regimen

## Give RATE
# ID1: two oral doses
# ID2: two IV doses
# ID3: oral then IV dose
nm <- data.frame(
  ID = c(1, 1, 2, 2, 3, 3),
  EVID = 1,
  CMT = c(1, 1, 2, 2, 1, 2),
  AMT = 100,
  TIME = c(0, 24, 0, 24, 0, 24),
  RATE = c(0, 0, 100, 100, 0, 200),
  DV = 0
)

reg_nm <- nm_to_regimen(nm)

dat_nm_rate <- sim(
  seed = 604,
  ode = mod,
  parameters = par,
  omega = omega,
  covariates = list(
    "WT" = 70,
    "AGE" = 14,
    "CYP2C19unknown" = 0,
    "CYP2C19a1a2" = 0,
    "CYP2C19a2a2" = 0,
    "CYP2C19a2a3" = 0,
    "CYP2C19a1a3" = 0,
    "CYP2C19a3a3" = 0
  ),
  n_ind = 3,
  regimen = reg_nm,
  t_obs = seq(0, 36, 0.1)
)

dat_nm_rate |>
  filter(comp == "obs") |>
  mutate(id = factor(id)) |>
  ggplot(aes(x = t, y = y, group = id, color = id)) +
  geom_point() +
  geom_line()

## NO RATE --> instant infusion
# ID1: two oral doses
# ID2: two IV doses
# ID3: oral then IV dose
nm <- data.frame(
  ID = c(1, 1, 2, 2, 3, 3),
  EVID = 1,
  CMT = c(1, 1, 2, 2, 1, 2),
  AMT = 100,
  TIME = c(0, 24, 0, 24, 0, 24),
  DV = 0
)

reg_nm <- nm_to_regimen(nm)

dat_nm <- sim(
  seed = 604,
  ode = mod,
  parameters = par,
  omega = omega,
  covariates = list(
    "WT" = 70,
    "AGE" = 14,
    "CYP2C19unknown" = 0,
    "CYP2C19a1a2" = 0,
    "CYP2C19a2a2" = 0,
    "CYP2C19a2a3" = 0,
    "CYP2C19a1a3" = 0,
    "CYP2C19a3a3" = 0
  ),
  n_ind = 3,
  regimen = reg_nm,
  t_obs = seq(0, 36, 0.1)
)

dat_nm |>
  filter(comp == "obs") |>
  mutate(id = factor(id)) |>
  ggplot(aes(x = t, y = y, group = id, color = id)) +
  geom_point() +
  geom_line()

Here, we give CMT with and without RATE, though in any data set with mixed oral/infusion dosing, you will have RATE.

With RATE, the infusion is not instantaneous:

Without RATE, the infusion is instant:

[roninsightrx]

What if we have two doses: on a first dose RATE=0 (bolus loading) and on second dose RATE=1? Then all doses will be treated as bolus, I think. Would add that as a test case to check, and potentially update behavior?

[dominic-irx]

The two blocks below are for backwards compatibility. If you have CMT defined (which you really should), then the previous if statement works fine. If I were writing this function from scratch, compartment would be a required column.

[roninsightrx]

looks good, but one last question about edge case

[roninsightrx]

lgtm