nlmixr2 / rxode2

rxode2
https://nlmixr2.github.io/rxode2/
GNU General Public License v3.0
28 stars 7 forks source link

Feature Request: Emperical Steady state #713

Open mattfidler opened 4 months ago

mattfidler commented 4 months ago

See:

https://nmhelp.tingjieguo.com/empiricalss.htm

I'm thinking that it is actually the same as what rxode2 does currently. Perhaps NONMEM may use some of the steady state routines inside of certain ODEs.

If this is allowed, it will allow you to specify empirical steady state max values for arbitrary doses, which isn't supported.

kylebaron commented 4 months ago

This is what mrgsolve does too. I always thought it was a bit of a workaround because I didn't know another way to do it. But Empirical Steady State sounds a lot better.

mattfidler commented 4 months ago

The difference here is you can change some options with addl. Monolix doesn't even check that that system has reached steady state.

You could possibly use root finding to get steady state, though it quite similar to the idea of the "emperical steady state"

Here is an example with the lsodar which is the rountine that would need to be used for this type of steady state:

library(deSolve)
## =======================================================================
## Example 2:
##   using lsodar to estimate steady-state conditions
## =======================================================================

## Bacteria (Bac) are growing on a substrate (Sub)
model <- function(t, state, pars) {
with (as.list(c(state, pars)), {
             ##        substrate uptake             death     respiration
             dBact <-  gmax*eff*Sub/(Sub+ks)*Bact - dB*Bact - rB*Bact
             dSub  <- -gmax    *Sub/(Sub+ks)*Bact + dB*Bact            + input

             return(list(c(dBact,dSub)))
             })
}

## root is the condition where sum of |rates of change|
## is very small

rootfun <- function (t, state, pars) {
dstate <- unlist(model(t, state, pars)) # rate of change vector
return(sum(abs(dstate)) - 1e-10)
}

pars <- list(Bini = 0.1, Sini = 100, gmax = 0.5, eff = 0.5,
                  ks = 0.5, rB = 0.01, dB = 0.01, input = 0.1)

tout    <- c(0, 1e10)
state   <- c(Bact = pars$Bini, Sub = pars$Sini)
out     <- lsodar(state, tout, model, pars, rootfun = rootfun)
print(out)
#>       time     Bact          Sub
#> 1    0.000 0.100000 100.00000000
#> 2 2583.657 3.333333   0.04347826

Created on 2024-04-19 with reprex v2.0.2

To me, the difference is switching between a true "root finding" ODE solving routine vs coding up the root finding internally (which is what we do in both rxode2 and mrgsolve).

I haven't read your code recently, though I do remember there are some subtle differences in root finding in things like steady state of continuous infusions. It may be more "stable" to find a root that is good-enough hence the introduction of the empirical steady state.

Since lsodar is in fortran I don't really feel the need to do anything other than empirical solutions (instead of true steady-state where the tolerance is probably 0).

mattfidler commented 4 months ago

In fact, I may simply change the empirical steady state records to a steady state record with some sort of warning.

mattfidler commented 3 months ago

This will use etTrans to translate negative addl records to steady state with warning.