Error: the complex F is not supported by babelmixr2

[john-harrold]

I'm testing the models in the model library (nlmixr2lib) trying to convert them to different formats. I'm having trouble converting this one to monolix.

library(babelmixr2)
#> Loading required package: nlmixr2
#> Loading required package: nlmixr2data
library(rxode2)
#> rxode2 2.1.3.9000 using 5 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
rx_fun = 
function() {
    covariateData <- list(WT = "Body weight in kg")
    description <- "Dupilumab PK model (Kovalenko 2020)"
    reference <- "Kovalenko P, Davis JD, Li M, et al. Base and Covariate Population Pharmacokinetic Analyses of Dupilumab Using Phase 3 Data. Clinical Pharmacology in Drug Development. 2020;9(6):756-767. doi:10.1002/cpdd.780"
    units <- list(time = "day", dosing = "mg")
    ini({
        lvc <- 0.908258560176891
        label("central volume (L)")
        lke <- -2.92994453301599
        label("elimination rate (1/d)")
        lkcp <- -1.54646311327271
        label("central-to-peripheral rate (1/d)")
        Mpc <- 0.686
        label("ratio of kcp and kpc (kpc is peripheral to central rate with units of 1/d)")
        lka <- -1.36257783450257
        label("absorption rate (1/d)")
        lMTT <- -2.25379492882461
        label("mean transit time (d)")
        lVm <- 0.0676586484738149
        label("maximum target-mediated rate of elimination (mg/L/d)")
        Km <- fix(0.01)
        label("Michaelis-Menten constant (mg/L)")
        lfdepot <- -0.441610554744518
        label("Bioavailability (fraction)")
        e_wt_vc <- 0.711
        label("Exponent of weight on central volume (unitless)")
        CcpropSd <- c(0, 0.15)
        label("Proportional residual error (fraction)")
        CcaddSd <- fix(0, 0.03)
        label("Additive residual error (mg/L)")
        etalvc ~ 0.192
        etalke ~ 0.285
        etalka ~ 0.474
        etalvm ~ 0.236
        etamtt ~ 0.525
    })
    model({
        vc <- exp(lvc + etalvc) * (WT/75)^e_wt_vc
        ke <- exp(lke + etalke)
        kcp <- exp(lkcp)
        ka <- exp(lka + etalka)
        MTT <- exp(lMTT + etamtt)
        Vm <- exp(lVm + etalvm)
        kpc <- kcp/Mpc
        ktr <- (3 + 1)/MTT
        d/dt(depot) <- -ktr * depot
        d/dt(transit1) <- ktr * (depot - transit1)
        d/dt(transit2) <- ktr * (transit1 - transit2)
        d/dt(transit3) <- ktr * transit2 - ka * transit3
        d/dt(central) <- ka * transit3 - ke * central - kcp * 
            central + kpc * periph - central * (Vm/(Km + central/vc))
        d/dt(periph) <- kcp * central - kpc * periph
        f(depot) <- exp(lfdepot)
        Cc <- central/vc
        Cc ~ add(CcaddSd) + prop(CcpropSd)
    })
}

rx_obj = rxode2::rxode2(rx_fun)
export_name = "test"

library(tibble)
dataset = 
tribble(
~id,  ~low,  ~time,  ~high,  ~cmt,  ~amt,  ~rate,  ~ii,  ~addl,  ~evid,  ~ss,  ~dur,  ~WT,  ~dv,
  1,    NA,      0,     NA,  'Cc',    NA,     NA,   NA,     NA,      0,   NA,    NA,    1,    0,
  1,    NA,      1,     NA,  'Cc',    NA,     NA,   NA,     NA,      0,   NA,    NA,    1,    0,
  2,    NA,      0,     NA,  'Cc',    NA,     NA,   NA,     NA,      0,   NA,    NA,    1,    0,
  2,    NA,      1,     NA,  'Cc',    NA,     NA,   NA,     NA,      0,   NA,    NA,    1,    0)

res = nlmixr2::nlmixr2(rx_obj, dataset, "monolix", babelmixr2::monolixControl(modelName=export_name, runCommand=NA))
#> Error : the complex F is not supported by babelmixr2
#> Error: the complex F is not supported by babelmixr2

^{Created on 2024-08-20 with reprex v2.1.1}

[mattfidler]

That is because it isn't supported because f(depot) is an equation. More parsing could probably support this particular expression, though I am not likely to get to this soon.

[mattfidler]

The error could be a bit more clear, I suppose.

[john-harrold]

Is there a way I can detect it's presence in rx_obj above?

[mattfidler]

Nope.

[mattfidler]

At least not as the nlmixr2/rxode2 interface is currently written.

[john-harrold]

Can I add that to the props? I know this isn't the place to request it :)

[mattfidler]

That is where it would likely go, though I believe I wouldn't support this directly since it is too much of a edge case for properties, rather I would support what properties (aka initial condition, f, rate, dur, alag) are present in the model and for what state).

Then you could use modelExtract() to parse this out yourself. I would say that if this changes, it wouldn't work so this idea is fragile. I think you should base it on if the model translates instead. But that is just my thoughts right now.

There are likely cases where a model won't import from nonmem or export to nonmem, and the same is true for monolix.

[john-harrold]

I'm sorry, you would rather support f, but isn't this f?

[john-harrold]

Oh I see. Can you add what properties are supported? I can easily use that.

[mattfidler]

I will add what properties are defined by the model.

For example in this model you would have some way of figuring out that f(depot) is in the model.

[john-harrold]

Awesome. Thank you.

[john-harrold]

If I wanted to conditionally skip tests for Monolix conversion where f is causing failure would I look at the Property column from the props$cmtProp data frame for a value of f?

$cmtProp
  Compartment Property
1       depot        f