transit function does not work with nlmixr with multiple administrations

[NicolaMelillo]

Hello,

I noticed that there seem to be some issues with the transit function when using nlmixr in case the dataset has multiple doses. It seems that only the first administration is considered when simulating the model.

Here it is a toy example with this issue.

library(nlmixr)
library(RxODE)

## load dataset
df.nm <- theo_md

## define model
my.model.pk <- function(){
  ini({

    lV   <- 0.2
    lka  <- -2
    ln   <- 0.3
    lCL  <- 0.9
    lktr <- 3

    eta.lka  ~ 0.5
    eta.lCL  ~ 0.5
    eta.ln   ~ 0.5

    err.pk <- 0.1

  })

  model({

    V   <- exp(lV)
    CL  <- exp(lCL + eta.lCL)
    ka  <- exp(lka + eta.lka) 
    ktr <- exp(lktr)
    n   <- exp(ln + eta.ln)

    mtt <- (n+1)/ktr
    bio <- 1

    d/dt(depot)   <-  transit(n,mtt,bio) - ka * depot
    d/dt(central) <-  ka * depot - CL/V * central

    central_conc <- central / V
    central_conc ~ add(err.pk)

  })

}

## fit the model
fit.nm <- nlmixr(my.model.pk, df.nm, est="saem",
                 control = saemControl(seed=2022,nBurn=500,nEm=500, transitAbs = TRUE))

## individual simulation
plot(augPred(fit.nm))

Individual fitting results are the following.

When using a dataset with single doses (e.g., theo_sd) everything works. When removing the transit absorption everything works as well.

I tried to see if there were problems with RxODE by simulating two doses, but everything seems to work.

library(RxODE)

mod <- RxODE({
  ## Table 3 from Savic 2007
  cl = 17.2 # (L/hr)
  vc = 45.1 # L
  ka = 0.3 # 1/hr
  mtt = 0.37 # hr
  bio=1
  n = 30
  k = cl/vc
  ktr = (n+1)/mtt
  d/dt(depot) = transit(n,mtt,bio)-ka*depot
  d/dt(cen) = ka*depot-k*cen
})

et.dose <- tibble(
  EVID = c(1,1),
  TIME = c(0,10),
  AMT  = c(10,10)
)

t.grid = seq(0,40,length.out=200)
et.obs <- tibble(
  EVID = numeric(length(t.grid)),
  TIME = t.grid,
  AMT  = numeric(length(t.grid))
)

et <- rbind(et.dose,et.obs)

transit <- rxSolve(mod, et)

plot(transit, cen, ylab="Central Concentration")

Simulated profiles are as follows.

Thanks, best wishes,

Nicola

[mattfidler]

Thank you for noting this. Since this can be identified using ODEs, perhaps the pure ode approach will work for now? Have you tried it?

[NicolaMelillo]

Hello Matt,

Thank you for your quick reply! For pure ode approach you mean expressing the transit function as explained here?

If this is the case, I have tried, but the results are like the ones with the transit function. Here it is the code.

library(nlmixr)
library(RxODE)

## load dataset
df.nm <- theo_md

## define model
my.model.pk <- function(){
  ini({

    lV   <- 0.2
    lka  <- -2
    ln   <- 0.3
    lCL  <- 0.9
    lktr <- 3

    eta.lka  ~ 0.5
    eta.lCL  ~ 0.5
    eta.ln   ~ 0.5

    err.pk <- 0.1

  })

  model({

    V   <- exp(lV)
    CL  <- exp(lCL + eta.lCL)
    ka  <- exp(lka + eta.lka) 
    ktr <- exp(lktr)
    n   <- exp(ln + eta.ln)

    # mtt <- (n+1)/ktr
    bio <- 1

    d/dt(depot)   <-  exp(log(bio*podo)+log(ktr)+n*log(ktr*t)-ktr*t-lgammafn(n+1)) - ka * depot
    d/dt(central) <-  ka * depot - CL/V * central

    central_conc <- central / V
    central_conc ~ add(err.pk)

  })

}

## fit the model
fit.nm <- nlmixr(my.model.pk, df.nm, est="saem",
                 control = saemControl(seed=2022,nBurn=500,nEm=500, transitAbs = TRUE))

## individual simulation
plot(augPred(fit.nm))

Individual fittings are as follows.

Perhaps you mean to directly write the ODEs of a certain number of transit compartments? I tried as well that approach, but the problem I have requires a BSV term on the number of absorption compartments, for this reason the function transit was ideal.

Thanks,

Nicola

[mattfidler]

Something similar: This worked for me with the CRAN nlmixr

library(nlmixr)
library(RxODE)
#> RxODE 1.1.4 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

## load dataset
df.nm <- theo_md

## define model
my.model.pk <- function(){
  ini({

    lV   <- 0.2
    lka  <- -2
    ln   <- 0.3
    lCL  <- 0.9
    lktr <- 3

    eta.lka  ~ 0.5
    eta.lCL  ~ 0.5
    eta.ln   ~ 0.5

    err.pk <- 0.1

  })

  model({

    V   <- exp(lV)
    CL  <- exp(lCL + eta.lCL)
    ka  <- exp(lka + eta.lka) 
    ktr <- exp(lktr)
    n   <- exp(ln + eta.ln)

    # mtt <- (n+1)/ktr
    bio <- 1
    if (amt != 0) tdos <- t
    td <- t - tdos
    d/dt(depot)   <-  exp(log(bio*podo)+log(ktr)+n*log(ktr*td)-ktr*td-lgamma(n+1)) - ka * depot
    d/dt(central) <-  ka * depot - CL/V * central

    central_conc <- central / V
    central_conc ~ add(err.pk)

  })

}

## fit the model
fit.nm <- nlmixr(my.model.pk, df.nm, est="saem",
                 control = saemControl(print=0))
#> → generate SAEM model
#> ✔ done
#> Calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> Error : Error calculating covariance via linearization
#> Calculating residuals/tables
#> done
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod
#> = FALSE, : Covariance matrix non-positive definite, corrected by sqrtm(fim %*%
#> fim)
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : Linearization of FIM could not be used to calculate covariance.
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : NAs introduced by coercion

## individual simulation
plot(fit.nm)

#> Warning: Transformation introduced infinite values in continuous x-axis
#> Warning in self$trans$transform(x): NaNs produced
#> Warning: Transformation introduced infinite values in continuous y-axis
#> Warning: Removed 2 rows containing missing values (geom_point).

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

#> Warning: Transformation introduced infinite values in continuous x-axis

^{Created on 2022-02-25 by the reprex package (v2.0.0)}

[NicolaMelillo]

Hi Matt,

That is great! The solution you provided works for me as well. Thank you!

Nicola

[mattfidler]

Glad to hear it.

Thank you for reporting the issue.

Matt

[NicolaMelillo]

Hi Matt,

I was playing a bit more with the pure ode approach. I noticed that it seems to ignore the part relative to the transit compartments. It seems that if transitAbs is not specified, the default is FALSE. In the toy example, if transitAbs=FALSE, the dose seems to be administered directly to the depot compartment.

library(nlmixr)
library(RxODE)

## load dataset
df.nm <- theo_md

## define model
my.model.pk <- function(){
  ini({

    lV   <- 0.2
    lka  <- -2
    ln   <- 0.3
    lCL  <- 0.9
    lktr <- 3

    eta.lka  ~ 0.5
    eta.lCL  ~ 0.5
    eta.ln   ~ 0.5

    err.pk <- 0.1

  })

  model({

    V   <- exp(lV)
    CL  <- exp(lCL + eta.lCL)
    ka  <- exp(lka + eta.lka) 
    ktr <- exp(lktr)
    n   <- exp(ln + eta.ln)

    mtt <- (n+1)/ktr
    bio <- 1

    if (amt != 0) tdos <- t
    td <- t - tdos
    d/dt(depot)   <-  exp(log(bio*podo)+log(ktr)+n*log(ktr*td)-ktr*td-lgamma(n+1)) - ka * depot
    d/dt(central) <-  ka * depot - CL/V * central

    central_conc <- central / V
    central_conc ~ add(err.pk)

  })

}

## fit the model
fit.nm <- nlmixr(my.model.pk, df.nm, est="saem",
                 control = saemControl(seed=2022,nBurn=200,nEm=200,transitAbs=FALSE))
plot(fit.nm)

In fact, if we multiply for 0 the transit compartments part of the equation, I obtain the same results of the code above.

d/dt(depot)   <-  exp(log(bio*podo)+log(ktr)+n*log(ktr*td)-ktr*td-lgamma(n+1)) * 0 - ka * depot
d/dt(central) <-  ka * depot - CL/V * central

If I set the option transitAbs=TRUE I obtain the following results.

I was playing a bit with the if conditions, I tried to put a condition on EVID as follows and something seems to work, although I got some errors...

library(nlmixr)
library(RxODE)

## load dataset
df.nm <- theo_md

## define model
my.model.pk <- function(){
  ini({

    lV   <- 0.2
    lka  <- -2
    ln   <- 0.3
    lCL  <- 0.9
    lktr <- 3

    eta.lka  ~ 0.5
    eta.lCL  ~ 0.5
    eta.ln   ~ 0.5

    err.pk <- 0.1

  })

  model({

    V   <- exp(lV)
    CL  <- exp(lCL + eta.lCL)
    ka  <- exp(lka + eta.lka) 
    ktr <- exp(lktr)
    n   <- exp(ln + eta.ln)

    mtt <- (n+1)/ktr
    bio <- 1

    if (EVID != 0) {
      tdos <- t
    } else{
      tdos <- 0
    }
    td <- t - tdos
    d/dt(depot)   <-  exp(log(bio*podo)+log(ktr)+n*log(ktr*td)-ktr*td-lgamma(n+1)) - ka * depot
    d/dt(central) <-  ka * depot - CL/V * central

    central_conc <- central / V
    central_conc ~ add(err.pk)

  })

}

## fit the model
fit.nm <- nlmixr(my.model.pk, df.nm, est="saem",
                 control = saemControl(seed=2022,nBurn=200,nEm=200,transitAbs=TRUE))
plot(fit.nm)

#Needed Covariates:
#[1] "EVID"
#Calculating residuals/tables
#Error : Column names `EVID.cov` and `EVID.cov` must not be duplicated.
#Use .name_repair to specify repair.
#Caused by error in `stop_vctrs()`:
#! Names must be unique.
#x These names are duplicated:
# * "EVID.cov" at locations 25, 26, and 27.
#Error in .fitFun(.ret) : FOCEi problem not allocated
#This can happen when sympy<->nlmixr interaction is not working correctly.
#Error in (function (data, inits, PKpars, model = NULL, pred = NULL, err = NULL,  : 
# Could not fit data.
#Warning message:
#In (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE,  :
# Error calculating nlmixr object, return classic object

Despite the errors, individual fits seems fine.

If I multiply for 0 the transit compartment part of the equation as per above, concentration stays always equal to 0.

I continue playing, I'll update you if I find something.

Nicola

[NicolaMelillo]

Hello again,

The above solution didn't work. Basically we would need to have the time after dose tad in the model section. Likely this is not the most elegant solution, but given that SAEM supports time dependent covariates, I tried to define tad as a model covariate and use it within the multicompartmental absorption function.

It does appear to work with me.

Code and results are as follows.

library(nlmixr)
library(RxODE)
library(tidyverse)

## load dataset
df.nm <- theo_md %>%
  mutate(
    TLAST_COV = ifelse(EVID!=0,TIME,NA)
  ) %>%
  fill(TLAST_COV) %>%
  mutate(TAD = TIME - TLAST_COV)

## define model
my.model.pk <- function(){
  ini({

    lV   <- 0.2
    lka  <- -2
    ln   <- 0.1
    lCL  <- 0.9
    lktr <- -1

    eta.lka  ~ 0.5
    eta.lCL  ~ 0.5
    eta.ln   ~ 0.5

    err.pk <- 0.1

  })

  model({

    V   <- exp(lV)
    CL  <- exp(lCL + eta.lCL)
    ka  <- exp(lka + eta.lka) 
    ktr <- exp(lktr)
    n   <- exp(ln + eta.ln)

    mtt <- (n+1)/ktr
    bio <- 1
    COV2 <- podo

    td <- TAD
    d/dt(depot)   <-  exp(log(bio*podo)+log(ktr)+n*log(ktr*td)-ktr*td-lgamma(n+1)) - ka * depot
    d/dt(central) <-  ka * depot - CL/V * central

    central_conc <- central / V
    central_conc ~ add(err.pk)

  })

}

## fit the model
fit.nm <- nlmixr(my.model.pk, df.nm, est="saem",
                 control = saemControl(seed=2022,nBurn=200,nEm=200,transitAbs=TRUE, calcTables = TRUE),
                 table = tableControl(addDosing = FALSE),
                 subsetNonmem = FALSE)

plot(fit.nm)

Does this make sense?

Nicola

[mattfidler]

This does make sense to me. In the future I may add an evid for transitAbs=TRUE. In theory both solutions should give similar values. In the case where the TAD is calculated in the model the ODE integration routine uses the correct time after dose for the points between the observations. In the case where it is included in the dataset, it uses LOCF integration between the time points.

Can you tell me why you think the first solution is not correct (I'm guessing it didn't run with your data)?

[NicolaMelillo]

Hi Matt,

So, for first solution you mean the one with the condition on EVID? In that case, I think it is wrong because td it is equal to tad. Probably it is due to the else condition. If I remove it it does not work. To avoid the error in the table calculation, it was better to use another covariate for the condition, equal to EVID but with a different name.

library(nlmixr)
library(RxODE)
library(tidyverse)

## load dataset
df.nm <- theo_md %>%
  mutate(
    COV_EV = EVID
  )

## define model
my.model.pk <- function(){
  ini({

    lV   <- 0.2
    lka  <- -2
    ln   <- 0.3
    lCL  <- 0.9
    lktr <- 3

    eta.lka  ~ 0.5
    eta.lCL  ~ 0.5
    eta.ln   ~ 0.5

    err.pk <- 0.1

  })

  model({

    V   <- exp(lV)
    CL  <- exp(lCL + eta.lCL)
    ka  <- exp(lka + eta.lka) 
    ktr <- exp(lktr)
    n   <- exp(ln + eta.ln)

    mtt <- (n+1)/ktr
    bio <- 1

    if (COV_EV != 0) {
      tdos <- t
    }else{
      tdos <- 0
    }
    td <- t - tdos
    d/dt(depot)   <-  exp(log(bio*podo)+log(ktr)+n*log(ktr*td)-ktr*td-lgamma(n+1)) - ka * depot
    d/dt(central) <-  ka * depot - CL/V * central

    central_conc <- central / V
    central_conc ~ add(err.pk)

  })

}

## fit the model
fit.nm <- nlmixr(my.model.pk, df.nm, est="saem",
                 control = saemControl(seed=2022,nBurn=200,nEm=200,transitAbs=TRUE))
plot(fit.nm)

# plot tad and td
df.plot <- fit.nm %>%
  select(TIME, tad, td) %>%
  gather(var, value, tad:td)

ggplot(df.plot, aes(TIME, value, color = var)) + 
  geom_line(size=2) + 
  xlab("TIME") + ylab(NULL)

[mattfidler]

I will have to explore the td issue a bit more.

I was thinking a bit more and this will be fine as long as the transit part is gone. Whenever a transit dose is administered the tad and dose are reset, so there is some of the dose that may accumulate but they may be reset with each dose. This happens quickly enough with the transit comparment is a surrogate for a lag time.

So, I think it looks OK.

Matt

[mattfidler]

In nlmixr2 this is being handled differently as a separate evid (evid=7) so you can use it in multiple modalities.

[NicolaMelillo]

Hi Matt,

Ok, thanks for the answer! My aim was to model a crossover study and the drug should be completely washed out from one dose to the other. So, my understanding is that, in that case, the solution with the TAD passed as a time-varying covariate should be fine.

What is not fine is when the second dose is administered while the first dose is still being absorbed, correct?

Looking forward to try nlmixr2! Thanks

[mattfidler]

I agree!

[NicolaMelillo]

Great, thank you very much Matt.