Receiving drastically different results when comparing ODEs with linCmt()

[WouterAh]

Hi @mattfidler, I am encountering an issue when comparing a 2-compartment model specified using ODEs with a 2-compartment model specified using linCmt(). When fitting these two models to the same data they produce different results. Below I pasted a simplified version of my code/models, using the simulated Infusion_2CPT data (my original data also only includes infusion data and is structured in a similar manner). Running this code reproduces the issue I am experiencing with my original model. This got me wondering whether one of the models presented here is not specified correctly in nlmixr syntax and brings me to the question of why switching from ODEs to linCmt() might in this case produce different results?

library(pacman)
pacman::p_load("nlmixr2", "nlmixr2est", "rxode2", "nlmixr2plot", "nlmixr2data", 
               "lotri", "nlmixr2extra", "xpose.nlmixr2")

data <- Infusion_2CPT

two_compartment_ODE <- function(){ 
  ini({ 
    lV1 <- 4.2 #log initial volume central compartment
    lCL <- 1.3 #log initial clearance central compartment
    lV2 <- 1.6 #log initial volume Peripheral compartment
    lQ <- 3.9 #log initial exchange rate between central and peripheral
    eta.V1 ~ 0.1 #IIV for V1
    eta.CL ~ 0.1 #IIV for CL
    eta.V2 ~ 0.1 #IIV for V2
    eta.Q ~ 0.1 #IIV for Q
    prop.err <- 0.2
    add.err <- 1
  })
  model({
    cl <- exp(lCL + eta.CL)
    q <- exp(lQ + eta.Q)
    v <- exp(lV1 + eta.V1)
    vp2 <- exp(lV2 + eta.V2)
    d/dt(cent) = -cl/v * cent - q/v * cent + q/vp2 * peri
    d/dt(peri) = q/v * cent - q/vp2 * peri
    cp = (cent/v)
    cp ~ add(add.err) + prop(prop.err)
  })
}

two_compartment_solved <- function(){ 
  ini({ 
    lV1 <- 4.2 #log initial volume central compartment
    lCL <- 1.3 #log initial clearance central compartment
    lV2 <- 1.6 #log initial volume Peripheral compartment
    lQ <- 3.9 #log initial exchange rate between central and peripheral
    eta.V1 ~ 0.1 #IIV for V1
    eta.CL ~ 0.1 #IIV for CL
    eta.V2 ~ 0.1 #IIV for V2
    eta.Q ~ 0.1 #IIV for Q
    prop.err <- 0.2
    add.err <- 1
  })
  model({
    cl <- exp(lCL + eta.CL)
    q <- exp(lQ + eta.Q)
    v <- exp(lV1 + eta.V1)
    vp2 <- exp(lV2 + eta.V2)
    cp = linCmt()
    cp ~ add(add.err) + prop(prop.err)
  })
}

fit.2cpt.ODE <- nlmixr2(two_compartment_ODE, data = data, est = "focei",
                        table=tableControl(cwres=T, npde=T))

fit.2cpt.solved <- nlmixr2(two_compartment_solved, data = data, est = "focei",
                           table=tableControl(cwres=T, npde=T))


print(fit.2cpt.ODE)
print(fit.2cpt.solved)

Best regards,
Wouter

[mattfidler]

Hi Wouter,

The linCmt() solutions have been verified by a variety of tests and these simulations have been reviewed and when run. While they produce different results, I don't believe they are drastically different.

These two tests were part of the original tests for the linCmt() solutions of focei. We saw differences but we didn't believe they were drastically different.

I am re-running these again to see if anything changed in the solutions.

As a note, the ODE solutions and linCmt() solutions will give slight different results because of their methods. This leads to different solutions when run. This is the same as many other nonlinear mixed effects modeling tools.

[mattfidler]

Looking at the models above again, we ran the proportional errors only.

[WouterAh]

Hi Matthew,

Thanks for your very fast response! I am indeed aware that slight differences in results are to be expected when comparing ODEs with linCmt(), however, when I was comparing these two types of model specification some parameter estimates turned out to increase/decrease by more than 25%. I appreciate you re-running these tests to see if you can reproduce the differences I find between specifying a two-compartment model with ODEs and linCmt().

Best,
Wouter

(p.s. I see now that I made a mistake in my initial parameter values and meant to switch around initial values for V2 and Q, doing so makes the population parameter estimates a lot more comparable but predictions and individual parameter estimates still seem to be way off for the solved model)

[mattfidler]

I can see there are differences between the two, but I'm unsure why. We have many tests that test for equality of the solved and ODE systems.

The only thing I can think of is that perhaps using a multi-exponential rather than an advan solution would possibly do a bit better.

This is only a possibility and takes alot of time to code, so I'm not going to immediately do anything for except for file a bug in rxode2 for these models.

[WouterAh]

Thanks again, I understand that you won't immediately be able to do anything about this issue. For now I'll continue using ODEs for my model development and I'll keep an eye out on the RxODE github to follow any progress.

[mattfidler]

I see another advan style solution. Perhaps it reduces calculation time and complexity. If you are tracking the issue it is here.

nlmixr2/rxode2#441

[WouterAh]

Thanks for coming back to me, I haven't had the time to read through it in depth but it certainly looks interesting as I also want to incorporate a non-zero initial condition in the model that I am working on. I'll make sure to follow that issue on the rxode2 GitHub as well and am very curious to see if this alternative analytical solution can at some point be incorporated in nlmixr2/rxode2!

[mattfidler]

In theory the non-zero initial condition currently works, if the system is stable. The current solution matches NONMEM's output (based on nonmem2rx check so far).