Error - Drug/metabolite modeing after zero oder absorption.

[Lee21A]

I tried to fit Parent and metabolite simultaneously after 0 order absorption of parent dug based on 1-compartment model.

km: conversion rate from parent to metabolite
cl : clearance for parent
clm : clearance for metabolite
v : volume of distribution
dur0: duration time of 0 order absorption

Step 1. I made a model like this.

library(nlmixr2)
one.0abs.PM <- function() {
  ini({
    tkm <- log(c(0,1)) 
    tcl <- log(c(0,2)) 
    tclm <- log(c(0,10))  
    tv <- log(c(0,10))  
    tdur0 <- log(c(0,1))  
    eta.cl ~ 1
    eta.v ~ 1
    eta.km ~ 1 
    eta.clm ~ 1
    eta.dur0 ~ 1
    prop.err <- 1 
    prop.err2 <- 1 
  })
  
  model({
    km <- exp(tkm + eta.km)
    cl <- exp(tcl + eta.cl)
    clm <- exp(tclm + eta.clm)
    v <- exp(tv + eta.v)
    dur0    <- exp(tdur0+eta.dur0)
    d/dt(center) = - cl / v * center - km * center
    dur(center)=dur0
    d/dt(meta) = km * center - clm/ v * meta
    center(0)= 0
    meta(0)= 0
    cp = center / v
    cm = meta / v
    cp ~ prop(prop.err)   
    cm ~  prop(prop.err2) 
  })
}
nlmixr2(one.0abs.PM)

The result was fine, and information for multiple endpoint was

── Multiple Endpoint Model ($multipleEndpoint): ──  
  variable               cmt               dvid*
1  cp ~ … cmt='cp' or cmt=3 dvid='cp' or dvid=1
2  cm ~ … cmt='cm' or cmt=4 dvid='cm' or dvid=2
  * If dvids are outside this range, all dvids are re-numered sequentially, ie 1,7, 10 becomes 1,2,3 etc

Step 2. I made a data file considering the upper information.

I used numbers for CMT and DVID based on the upper multiple endpoint model.

Step 3. I tried to fit.

dat <- read.csv("0ABS_PM_Report3.csv",header = T,encoding = "UTF-8", na.strings=c("."))
fit <- nlmixr2(one.0abs.PM, dat, list(print=0), est="saem")

Step 4. The fitting result was weird, and there was a wrong assignment for DV data. (for example 3.8)

For every ID, DV (such as 3.8) of time 0.5 assigned as DVID=1 is “cm”, not “cp”.
This weird phenomenon applied for only time 0.5.

Step 5. I changed the model like this and the CMT info, But it was the same result.

cp ~ prop(prop.err)   |center  --> changed the number of CMT 3 to 1
cm ~  prop(prop.err2) |meta -->  changed the number of CMT 4 to 2

Step 5. I tried to change the info of DVID of CMT to “cp”, “cm” or “center”, “meta”. But the result was the same or I encountered an error like this.

'DVID'/'CMT' translation:
'CMT': 4
Error : 'dvid'->'cmt' or 'cmt' on observation record on a undefined compartment (use 'cmt()' 'dvid()') id: 1 row: 0
Error: 'dvid'->'cmt' or 'cmt' on observation record on a undefined compartment (use 'cmt()' 'dvid()') id: 1 row: 0

Could somebody help me?
I attached csv files I used.
Thanks in advance.

Lee
0ABS_PM_Report3.csv

[mattfidler]

This seems to happen in some edge case where the values for different end-points are tied and they are not synced appropriately. nlmixr2/rxode2#581 has a more precise reproducible example that shows this. While all the times are tied for some reason it only applies to the 0.5 time point.

In the mean-time, by simply offsetting the TIME by a very small amount for one of the compartments will make the fit perform a bit better, and the time-point at approximately 0.5 shows the right compartment:

library(nlmixr2est)
#> Warning: package 'nlmixr2est' was built under R version 4.3.1
#> Loading required package: nlmixr2data
#> Warning: package 'nlmixr2data' was built under R version 4.3.1

one.0abs.PM <- function() {
  ini({
    tkm <- log(c(0,1)) 
    tcl <- log(c(0,2)) 
    tclm <- log(c(0,10))  
    tv <- log(c(0,10))  
    tdur0 <- log(c(0,1))  
    eta.cl ~ 1
    eta.v ~ 1
    eta.km ~ 1 
    eta.clm ~ 1
    eta.dur0 ~ 1
    prop.err <- 1 
    prop.err2 <- 1 
  })
  
  model({
    km <- exp(tkm + eta.km)
    cl <- exp(tcl + eta.cl)
    clm <- exp(tclm + eta.clm)
    v <- exp(tv + eta.v)
    dur0    <- exp(tdur0+eta.dur0)
    d/dt(center) = - cl / v * center - km * center
    dur(center)=dur0
    d/dt(meta) = km * center - clm/ v * meta
    center(0)= 0
    meta(0)= 0
    cp = center / v
    cm = meta / v
    cp ~ prop(prop.err)   
    cm ~  prop(prop.err2) 
  })
}

dat <- read.csv("c:/tmp/0ABS_PM_Report3.csv",header = T,encoding = "UTF-8", na.strings=c("."))
dat$TIME[dat$CMT==4] <- dat$TIME[dat$CMT==4]+0.0001

fit <- nlmixr2(one.0abs.PM, dat, list(print=0), est="saem")
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> rxode2 2.0.13 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> Calculating covariance matrix
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem predOnly model 0...
#> → optimizing duplicate expressions in saem predOnly model 0...
#> → finding duplicate expressions in saem predOnly model 1...
#> → optimizing duplicate expressions in saem predOnly model 1...
#> → finding duplicate expressions in saem predOnly model 2...
#> → optimizing duplicate expressions in saem predOnly model 2...
#> ✔ done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 4128
#> → compress phiM in nlmixr2 object, save 77096
#> → compress parHist in nlmixr2 object, save 20704
#> → compress saem0 in nlmixr2 object, save 38696

print(fit)
#> ── nlmixr² SAEM OBJF by FOCEi approximation ──
#> 
#>  Gaussian/Laplacian Likelihoods: AIC() or $objf etc. 
#>  FOCEi CWRES & Likelihoods: addCwres() 
#> 
#> ── Time (sec $time): ──
#> 
#>         setup saem table compress other
#> elapsed 0.002 6.94  0.03     0.03 4.088
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#> 
#>              Est.    SE %RSE Back-transformed(95%CI)  BSV(CV%) Shrink(SD)%
#> tkm       -0.0368 0.224  611      0.964 (0.621, 1.5)   0.00970      94.6% 
#> tcl         -1.34 0.967 72.2    0.262 (0.0393, 1.74)     0.124      95.9% 
#> tclm       0.0805 0.236  294      1.08 (0.682, 1.72)   0.00370      94.9% 
#> tv          0.163 0.188  115       1.18 (0.814, 1.7) 0.0000599      97.1% 
#> tdur0     -0.0984 0.144  147      0.906 (0.683, 1.2)    0.0335      98.4% 
#> prop.err    0.594                              0.594                      
#> prop.err2   0.332                              0.332                      
#>  
#>   Covariance Type ($covMethod): linFim
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Information about run found ($runInfo):
#>    • Undefined 'dvid' integer values in data: 1, 2 
#>    • package 'rxode2' was built under R version 4.3.1 
#>   Censoring ($censInformation): No censoring
#> 
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 75 × 25
#>   ID      TIME CMT      DV         PRED      RES   IPRED     IRES  IWRES  eta.cl
#>   <fct>  <dbl> <fct> <dbl>        <dbl>    <dbl>   <dbl>    <dbl>  <dbl>   <dbl>
#> 1 1     0.0001 cm      0   0.0000000451 -4.51e-8 4.51e-8 -4.51e-8 -3.01  5.08e-5
#> 2 1     0.5    cp      3.8 3.54          2.65e-1 3.54e+0  2.65e-1  0.126 5.08e-5
#> 3 1     0.500  cm      0.9 0.802         9.82e-2 8.02e-1  9.82e-2  0.369 5.08e-5
#> # ℹ 72 more rows
#> # ℹ 15 more variables: eta.v <dbl>, eta.km <dbl>, eta.clm <dbl>,
#> #   eta.dur0 <dbl>, cp <dbl>, cm <dbl>, center <dbl>, meta <dbl>, km <dbl>,
#> #   cl <dbl>, clm <dbl>, v <dbl>, dur0 <dbl>, tad <dbl>, dosenum <dbl>

print(fit %>% as.data.frame() %>% dplyr::mutate(TIME= round(TIME,3)) %>%
        dplyr::filter(TIME==0.5) %>%
        dplyr::select(ID, TIME, CMT, IPRED, DV))
#>   ID TIME CMT     IPRED  DV
#> 1  1  0.5  cp 3.5352614 3.8
#> 2  1  0.5  cm 0.8017849 0.9
#> 3  2  0.5  cp 3.5352739 3.4
#> 4  2  0.5  cm 0.8017826 0.7
#> 5  3  0.5  cp 3.5352973 2.8
#> 6  3  0.5  cm 0.8017842 1.1

^{Created on 2023-09-11 with reprex v2.0.2}

While it isn't an ideal fix, it allows you to use the CRAN version of nlmixr2 and still get reasonable results.

I am hoping to have a CRAN release of the nlmixr2 family of packages soon once a few things get fixed (including this particular issue now).

Thank you for raising this issue so we can fix this case of 2 tied endpoints not synced appropriately.

[mattfidler]

Since this is no longer the issue (and I have a different issue tagged for this problem). I will close this for now and tag you in the rxode2 issue so you can track progress.

[mattfidler]

@Lee21A this has been updated in development.

For the ODE system you are using it should work just fine.

This means with the new r-universe build you should be able to run without issue.

This would update everything needed; You will likely need to wait for an hour or so before the build is complete.

install.packages(c("dparser", "nlmixr2data", "lotri", "rxode2ll",
                   "rxode2parse", "rxode2random", "rxode2et",
                   "rxode2", "nlmixr2est", "nlmixr2extra", "nlmixr2plot",
                   "nlmixr2"),
                 repos = c('https://nlmixr2.r-universe.dev',
                           'https://cloud.r-project.org'))

[Lee21A]

I really appreciate your quick answer. It works!

[mattfidler]

Great!