Help please: dual absorption nlmixr model (first+zero)

[yjung1]

Hello, I was trying to fit a dual absorption model using the dataset that shows dual absorption. I have tried:

1. Zero-order absorption followed by the first-order absorption
2. First-order absorption followed by zero-order absorption

The first model (Zero + first) seems to work fine, but somehow, the second model does not recognize parameter ALAG2 So, I also made another dataset that includes both ALAG1 and ALAG2 and fitted to the first model (Zero+first) and it seems to estimate the ALAG2 parameter.

I have been stuck with this problem forever. Would be so appreciated if I could get some help, please. I have attached the dataset as well.

This is the first-model;

#### 1-comp Dual (Zero+first) absorption and First-order elimination ####

dat <- read.csv("1-comp_zf.csv") %>% 
  filter(time==0) %>% 
  mutate(AMT = 100000, RATE = -2, CMT = 2) 

data <- read.csv("1-comp_zf.csv") %>% 
  filter(time==0) %>% 
  mutate(AMT = 100000, RATE = 0, CMT = 1) %>% 
  rbind(dat)

data_dual <- read.csv("1-comp_zf.csv") %>% 
  mutate(AMT = 0, RATE = 0, CMT = 2) %>% 
  rbind(data) %>% 
  filter(time %in% c(0, 1, 2, 4, 6, 8, 12, 16, 24)) %>% 
  rename(DV = CP, TIME = time) %>% 
  arrange(ID,TIME,RATE,CMT)

Z301 <- function() {
  ini({
    tka    <- log(c(0,1))
    tcl    <- log(c(0,10))
    tv     <- log(c(0,100))
    tD2    <- log(c(0,2))
    tf1    <- log(c(0,0.5,1))
    talag1 <- log(c(0,6))

    eta.ka ~ 0.1
    eta.cl ~ 0.1
    eta.v  ~ 0.1

    prop.err <- 0.1
  })

  model({
    KA    <- exp(tka + eta.ka)
    CL    <- exp(tcl + eta.cl)
    VC    <- exp(tv  + eta.v)
    D2    <- exp(tD2)
    F1    <- exp(tf1)
    ALAG1 <- exp(talag1)

    K20 = CL/VC
    S2 = VC

    d/dt(A1) = -KA*A1
    alag(A1) = ALAG1
    d/dt(A2) = KA*A1 - K20*A2
    dur(A2)  = D2
    f(A1) <- F1
    f(A2) <- 1-F1
    CP = A2/S2
    CP ~ prop(prop.err)

  })
}

Second model I am stuck with...

#### 1-comp Dual (First+zero) absorption and First-order elimination ####

dat <- read.csv("1-comp_zf.csv") %>% 
  filter(time==0) %>% 
  mutate(AMT = 100000, RATE = -2, CMT = 2) 

data <- read.csv("1-comp_zf.csv") %>% 
  filter(time==0) %>% 
  mutate(AMT = 100000, RATE = 0, CMT = 1) %>% 
  rbind(dat)

data_dual <- read.csv("1-comp_zf.csv") %>% 
  mutate(AMT = 0, RATE = 0, CMT = 2) %>% 
  rbind(data) %>% 
  filter(time %in% c(0, 1, 2, 4, 6, 8, 12, 16, 24)) %>% 
  rename(DV = CP, TIME = time) %>% 
  arrange(ID,TIME,CMT, RATE)

Z201 <- function() {
  ini({
    tka    <- log(c(0,3))
    tcl    <- log(c(0,5))
    tv     <- log(c(0,100))
    tD2    <- log(c(0,1))
    tf1    <- log(c(0,0.5,1))
    talag2 <- log(c(0,6))

    eta.ka ~ 0.1
    eta.cl ~ 0.1
    eta.v  ~ 0.1

    prop.err <- 0.1
  })

  model({
    KA    <- exp(tka + eta.ka)
    CL    <- exp(tcl + eta.cl)
    VC    <- exp(tv + eta.v)
    D2    <- exp(tD2)
    F1    <- exp(tf1)
    ALAG2 <- exp(talag2)

    K20 = CL/VC
    S2 = VC

    d/dt(A1) = -KA*A1
    d/dt(A2) = KA*A1 - K20*A2
    dur(A2)  = D2
    alag(A2) = ALAG2
    f(A1) <- F1
    f(A2) <- 1-F1
    CP = A2/S2
    CP ~ prop(prop.err)

  })
}

dataset.zip

[yjung1]

The dataset I used was simulated using the parameters I have given temporarily; Dose = 100mg CL = 10 L/h V = 100 L Ka = 1 /h ALAG1 = 6 (Zero + first) ALAG2 = 5 (First + zero)

BSV on CL, V, Ka Residual error: Proportional error

[mattfidler]

It is very odd. I looked at the underlying RxODE and the results between the two models are are the same both in terms of translation and solving.

library(RxODE)
#> RxODE 1.1.0 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
library(tidyverse)
library(testthat)
#> 
#> Attaching package: 'testthat'
#> The following object is masked from 'package:dplyr':
#> 
#>     matches
#> The following object is masked from 'package:purrr':
#> 
#>     is_null
#> The following object is masked from 'package:tidyr':
#> 
#>     matches

setwd("~/src/issues")

#### 1-comp Dual (Zero+first) absorption and First-order elimination ####

dat <- read.csv("1-comp_zf.csv") %>%
  filter(time == 0) %>%
  mutate(AMT = 100000, RATE = -2, CMT = 2)

data <- read.csv("1-comp_zf.csv") %>%
  filter(time == 0) %>%
  mutate(AMT = 100000, RATE = 0, CMT = 1) %>%
  rbind(dat)

data_dual <- read.csv("1-comp_zf.csv") %>%
  mutate(AMT = 0, RATE = 0, CMT = 2) %>%
  rbind(data) %>%
  filter(time %in% c(0, 1, 2, 4, 6, 8, 12, 16, 24)) %>%
  rename(DV = CP, TIME = time) %>%
  arrange(ID, TIME, RATE, CMT)

rx <- RxODE({
  param(tka, tcl, tv, tD2, tf1, talag1)
  cmt(A1)
  cmt(A2)
  KA ~ exp(tka)
  CL ~ exp(tcl)
  VC ~ exp(tv)
  D2 ~ exp(tD2)
  F1 ~ exp(tf1)
  ALAG1 ~ exp(talag1)
  K20 ~ CL / VC
  S2 ~ VC
  d / dt(A1) <- -KA * A1
  alag(A1) <- ALAG1
  d / dt(A2) <- KA * A1 - K20 * A2
  dur(A2) <- D2
  f(A1) <- F1
  f(A2) <- 1 - F1
  CP ~ A2 / S2
  nlmixr_pred <- CP
  cmt(CP)
})

data_dual2 <- read.csv("1-comp_zf.csv") %>%
  mutate(AMT = 0, RATE = 0, CMT = 2) %>%
  rbind(data) %>%
  filter(time %in% c(0, 1, 2, 4, 6, 8, 12, 16, 24)) %>%
  rename(DV = CP, TIME = time) %>%
  arrange(ID, TIME, CMT, RATE)

expect_equal(
  etTrans(data_dual, rx),
  etTrans(data_dual2, rx)
)

expect_equal(
  rxSolve(rx, data_dual, c(
    tka = -0.474897028364044, tcl = 2.32960078767426, tv = 4.56694058830605,
    tD2 = 0.424343867676095, tf1 = -0.896040171664832, talag1 = 1.90071932601878,
    prop.err = 0.093016424799712
  ), addDosing = TRUE, returnType = "data.frame"),
  rxSolve(rx, data_dual2, c(
    tka = -0.474897028364044, tcl = 2.32960078767426, tv = 4.56694058830605,
    tD2 = 0.424343867676095, tf1 = -0.896040171664832, talag1 = 1.90071932601878,
    prop.err = 0.093016424799712
  ), addDosing = TRUE, returnType = "data.frame")
)

^{Created on 2021-06-14 by the reprex package (v2.0.0)}

[mattfidler]

Running both models with focei gives the same results.

[mattfidler]

So, I'm assuming you are talking about the saem algorithm, right?

[mattfidler]

I didn't immediatly see the difference. The difference is in the ETA ka

[mattfidler]

Which algorithm is there a problem with so I can narrow down my debuging.

[yjung1]

the dataset is the same for both models, but the first one is the zero-order + first-order absorption which requires ALAG1, whereas model 2 (first-order + zero-order absorption) needing ALAG2 parameter. I'm not getting the same results, or model 2 does not seem to fit at all. I used the focei algorithm

[mattfidler]

model 2 doesn't work well with saem either.

[yjung1]

So does that mean nlmixr cannot build the first-order + zero-order (dual) absorption structure?

[yjung1]

I'm trying to figure out what I did wrong... are the codes wrong or do you think it's something to do with zero-order absorption and lag time not being recognized? Another topic but might be related? 1-compartment, zero-order absorption with lag time structure I'm getting initial estimates as the final estimates. So if I run the model with an initial of (0,4) for D1, I get 4 as the final estimate, an initial of (0,5) gives me 5 as the final estimate, and so on....

[mattfidler]

In model 2, you are trying to estimate the start and stop of the infusion combined with a first order absorption starting at time zero.

In model 1, you are trying to estimate the length of the infusion and the start of the first order absorption.

I believe that in general model 2 is more difficult to estimate regardless of what method you use. When you try to estimate the start of the infusion as well as the end combined with the first order absorption split it is likely not identifiable regardless of if you simulated it from an ODE. When you estimate the start time of the infusion it is confounded by the end time of the infusion, the estimates for each are highly interrelated.

In general, the time based dosing changes like duration of infusion and lag time are more difficult to estimate. This is why you are also getting no information for the different lag times.

When I test the system with a simple lag time everything is fine in RxODE

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

#### 1-comp Dual (Zero+first) absorption and First-order elimination ####

rx <- RxODE({
  K20 = CL/VC
  S2 = VC
  d/dt(A1) = -KA*A1
  d/dt(A2) = KA*A1 - K20*A2
  dur(A2)  = D2
  alag(A2) = ALAG2
  CP = A2/S2
})

t <- c(KA = 3,
  CL = 5,
  VC=100,
  D2=1,
  ALAG2=1)

et <- et(amt = 100000, rate = -2, cmt = 2)

a <- rxSolve(rx, et, t)

plot(a, A1, A2, CP) %>% print()

^{Created on 2021-06-16 by the reprex package (v2.0.0)}

[mattfidler]

Maybe there is some interaction between the two that causes the strange behavior (or the system is unidentifiable)

[mattfidler]

There seems to be an error solving in RxODE when mixing all of these together.

I have a simpler reproducible example below:

library(RxODE)
#> RxODE 1.1.0 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
library(ggplot2)
library(patchwork)

rx1 <- RxODE({
  K20 = CL/VC
  S2 = VC
  d/dt(A1) = -KA*A1
  d/dt(A2) = KA*A1 - K20*A2
  dur(A2)  = D2
  alag(A2) = ALAG2
  CP = A2/S2
})

rx2 <- RxODE({
  K20 = CL/VC
  S2 = VC
  d/dt(A1) = -KA*A1
  d/dt(A2) = KA*A1 - K20*A2
  dur(A2)  = D2
  alag(A2) = ALAG2
  f(A1) <- F1
  f(A2) <- 1 - F1
  CP = A2/S2
})

t <- c(KA = 3,
  CL = 5,
  VC=100,
  D2=4,
  F1=0,
  ALAG2=4)

et1 <- et(amt = 100000, rate = -2, cmt = 2)
et2 <- et1 %>%
  et(amt = 100000, cmt = 1)

a <- rxSolve(rx1, et1, t)
b <- rxSolve(rx2, et2, t)

g1 <- plot(a, A1, A2, CP) + ggtitle("direct")
g2 <- plot(b, A1, A2, CP) + ggtitle("mixed")

g1 / g2

^{Created on 2021-06-16 by the reprex package (v2.0.0)}

[mattfidler]

I did fix the issue in RxODE described above. A simple rerun didn't create a better fit, though.

You could try reinstalling RxODE and nlmixr from github and then trying with your data:

.remotes::install_github("nlmixrdevelopment/RxODE")
remotes::install_github("nlmixrdevelopment/nlmixr")

I haven't had time to see if I can pursue the issue further, though. I will be unable to look into it further until next week.

Here is the correct behavior based on the development version of RxODE

library(RxODE)
#> RxODE 1.1.0 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
library(ggplot2)
library(patchwork)

rx1 <- RxODE({
  K20 = CL/VC
  S2 = VC
  d/dt(A1) = -KA*A1
  d/dt(A2) = KA*A1 - K20*A2
  dur(A2)  = D2
  alag(A2) = ALAG2
  CP = A2/S2
})

rx2 <- RxODE({
  K20 = CL/VC
  S2 = VC
  d/dt(A1) = -KA*A1
  d/dt(A2) = KA*A1 - K20*A2
  dur(A2)  = D2
  alag(A2) = ALAG2
  f(A1) <- F1
  f(A2) <- 1 - F1
  CP = A2/S2
})

t <- c(KA = 3,
  CL = 5,
  VC=100,
  D2=4,
  F1=0,
  ALAG2=4)

et1 <- et(amt = 100000, rate = -2, cmt = 2)
et2 <- et1 %>%
  et(amt = 100000, cmt = 1)

a <- rxSolve(rx1, et1, t)
b <- rxSolve(rx2, et2, t)

g1 <- plot(a, A1, A2, CP) + ggtitle("direct")
g2 <- plot(b, A1, A2, CP) + ggtitle("mixed")

g1 / g2

^{Created on 2021-06-21 by the reprex package (v2.0.0)}

[yjung1]

thank you for your help. I will try re-installing!

[liangqiongyue]

1. Zero-order absorption followed by the first-order absorption

hi, for the first model (zero+first), why isn't the code written this way?

model({
    KA    <- exp(tka + eta.ka)
    CL    <- exp(tcl + eta.cl)
    VC    <- exp(tv  + eta.v)
    D2    <- exp(tD2)
    F1    <- exp(tf1)
    ALAG1 <- D2

    K20 = CL/VC
    S2 = VC

    d/dt(A1) = -KA*A1
    alag(A1) = ALAG1
    d/dt(A2) = KA*A1 - K20*A2
    dur(A2)  = D2
    f(A1) <- F1
    f(A2) <- 1-F1
    CP = A2/S2
    CP ~ prop(prop.err)

  })

That is to change the code to ALAG1 <- D2

Because don't you have to wait for zero-order absorption to end before you start first-order absorption?

[mattfidler]

Because don't you have to wait for zero-order absorption to end before you start first-order absorption?

There is no requirement that the zero-order absorption ends before the start of a first order absorption.