Requesting an example converting rxode2 to NONMEM

[john-harrold]

I've looked through the documentation for rxToNonmem() and I'm not sure what the first input is :).

Say I have the model below. How do I convert that into NONMEM or Monolix format? I'm assuming the
the ui argument is the rx below.

mod1 <- function() {
  ini(https://nlmixr2.github.io/rxode2/reference/ini.html)({
    # central 
    KA=2.94E-01
    CL=1.86E+01
    V2=4.02E+01
    # peripheral
    Q=1.05E+01
    V3=2.97E+02
    # effects
    Kin=1
    Kout=1
    EC50=200 
  })
  model(https://nlmixr2.github.io/rxode2/reference/model.html)({
    C2 <- centr/V2
    C3 <- peri/V3
    d/[dt](https://rdrr.io/r/stats/TDist.html)(depot) <- -KA*depot
    d/[dt](https://rdrr.io/r/stats/TDist.html)(centr) <- KA*depot - CL*C2 - Q*C2 + Q*C3
    d/[dt](https://rdrr.io/r/stats/TDist.html)(peri)  <- Q*C2 - Q*C3
    eff(0) <- 1
    d/[dt](https://rdrr.io/r/stats/TDist.html)(eff)   <- Kin - Kout*(1-C2/(EC50+C2))*eff
  })
}

rx = rxode2(mod1)

[mattfidler]

The way to create (and run) NONMEM from babelmixr2 is:

fit <- nlmixr2(rx, data, "nonmem")

You have to have a fully specified model with endpoints for the conversion to work.

Since the data is translated to NONMEM compatible data, and the estimation method is written in a specific way, there is no way to translate an arbitrary rxode2 code to a NONMEM control stream.

The function you found is really to convert parts of rxode2 to NONMEM code.

However, there is a stub within rxode2 ui that you can use: $nonmeModel. This stub only works when babelmixr2 is loaded and will use the default options:

library(babelmixr2)
#> Loading required package: nlmixr2
#> Loading required package: nlmixr2data

one.compartment <- function() {
  ini({
    tka <- log(1.57); label("Ka")
    tcl <- log(2.72); label("Cl")
    tv <- log(31.5); label("V")
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.sd <- 0.7
  })
  # and a model block with the error specification and model specification
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    d/dt(depot) <- -ka * depot
    d/dt(center) <- ka * depot - cl / v * center
    cp <- center / v
    cp ~ add(add.sd)
  })
}

f <- one.compartment()

message(f$nonmemModel)
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> $PROBLEM  translated from babelmixr2
#> ; comments show mu referenced model in ui$getSplitMuModel
#> 
#> $DATA one.compartment.csv IGNORE=@
#> 
#> $INPUT ID TIME EVID AMT DV CMT RXROW
#> 
#> $SUBROUTINES ADVAN13 TOL=6 ATOL=12 SSTOL=6 SSATOL=12
#> 
#> $MODEL NCOMPARTMENTS=2
#>      COMP(DEPOT, DEFDOSE) ; depot
#>      COMP(CENTER) ; center
#> 
#> $PK
#>   MU_1=THETA(1)
#>   MU_2=THETA(2)
#>   MU_3=THETA(3)
#>   KA=DEXP(MU_1+ETA(1)) ; ka <- exp(tka)
#>   CL=DEXP(MU_2+ETA(2)) ; cl <- exp(tcl)
#>   V=DEXP(MU_3+ETA(3)) ; v <- exp(tv)
#> 
#> $DES
#>   DADT(1) = - KA*A(1) ; d/dt(depot) = -ka * depot
#>   DADT(2) = KA*A(1)-CL/V*A(2) ; d/dt(center) = ka * depot - cl/v * center
#>   CP=A(2)/V ; cp = center/v
#> 
#> $ERROR
#>   ;Redefine LHS in $DES by prefixing with on RXE_ for $ERROR
#>   RXE_CP=A(2)/V ; cp = center/v
#>   RX_PF1=RXE_CP ; rx_pf1 ~ cp
#>   ; Write out expressions for ipred and w
#>   RX_IP1 = RX_PF1
#>   RX_P1 = RX_IP1
#>   W1=DSQRT((THETA(4))**2) ; W1 ~ sqrt((add.sd)^2)
#>   IF (W1 .EQ. 0.0) W1 = 1
#>   IPRED = RX_IP1
#>   W     = W1
#>   Y     = IPRED + W*EPS(1)
#> 
#> $THETA (0.45108   ) ; 1 - tka   
#>        (1.0006    ) ; 2 - tcl   
#>        (3.45      ) ; 3 - tv    
#>        (0,     0.7) ; 4 - add.sd
#> 
#> $OMEGA 0.6 ; eta.ka
#> $OMEGA 0.3 ; eta.cl
#> $OMEGA 0.1 ; eta.v
#> 
#> $SIGMA 1 FIX
#> 
#> $ESTIMATION METHOD=1 INTER MAXEVALS=10000 SIGDIG=3 SIGL=12 PRINT=1 NOABORT
#> 
#> $COVARIANCE
#> 
#> $TABLE ID ETAS(1:LAST) OBJI FIRSTONLY ONEHEADER NOPRINT
#>      FORMAT=s1PE17.9 NOAPPEND FILE=one.compartment.eta
#> 
#> $TABLE ID TIME IPRED PRED RXROW ONEHEADER NOPRINT
#>     FORMAT=s1PE17.9 NOAPPEND FILE=one.compartment.pred

^{Created on 2024-05-24 with reprex v2.1.0}

[john-harrold]

That's cool. And I see that

f$monolixModel

Will give me a character object with the model in Monolix format.

I think that was exactly what I want.

I've got ubiquity creating models in rxode2 format and I just want to use babelmixr2 to also provide NONMEM and Monolix models as well.

[mattfidler]

Monolix is a bit more complex because it requires 2 files.

[mattfidler]

library(babelmixr2)
#> Loading required package: nlmixr2
#> Loading required package: nlmixr2data

one.compartment <- function() {
  ini({
    tka <- log(1.57); label("Ka")
    tcl <- log(2.72); label("Cl")
    tv <- log(31.5); label("V")
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.sd <- 0.7
  })
  # and a model block with the error specification and model specification
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    d/dt(depot) <- -ka * depot
    d/dt(center) <- ka * depot - cl / v * center
    cp <- center / v
    cp ~ add(add.sd)
  })
}

f <- one.compartment()

message(f$monolixModel)
#> DESCRIPTION:
#> model translated from `babelmixr2` and `nlmixr2` function one.compartment to one.compartment-monolix.txt
#> 
#> [LONGITUDINAL]
#> input={ka,cl,v}
#> 
#> PK:
#> ; Define compartments with administrations
#> compartment(cmt=1, amount=depot)
#> ; Define PK macros
#> depot(type=1, target=depot, Tlag=0, p=1)
#> 
#> EQUATION:
#>    ddt_depot = - ka*depot
#>    ddt_center = ka*depot-cl/v*center
#>    cp = center/v
#>    rx_pred_cp = cp
#> 
#> OUTPUT:
#> output={rx_pred_cp}

message(f$mlxtran)
#> <DATAFILE>
#> 
#> [FILEINFO]
#> file='one.compartment-monolix.csv'
#> delimiter = comma
#> header = {ID, TIME, EVID, AMT, DV, ADM, nlmixrRowNums}
#> 
#> [CONTENT]
#> ID = {use=identifier}
#> TIME = {use=time}
#> EVID = {use=eventidentifier}
#> AMT = {use=amount}
#> DV = {use=observation, name=rx_prd_cp, type=continuous}
#> ADM = {use=administration}
#> 
#> <MODEL>
#> 
#> [INDIVIDUAL]
#> input={ka_pop, omega_ka, cl_pop, omega_cl, v_pop, omega_v}
#> 
#> DEFINITION:
#> ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
#> cl = {distribution=logNormal, typical=cl_pop, sd=omega_cl}
#> v = {distribution=logNormal, typical=v_pop, sd=omega_v}
#> 
#> [LONGITUDINAL]
#> input={add__sd}
#> file='one.compartment-monolix.txt'
#> 
#> DEFINITION:
#> rx_prd_cp={distribution = normal, prediction = rx_pred_cp, errorModel=constant(add__sd)}
#> 
#> <FIT>
#> data={rx_prd_cp}
#> model={rx_prd_cp}
#> 
#> <PARAMETER>
#> ka_pop={value=1.57, method=MLE}
#> cl_pop={value=2.72, method=MLE}
#> v_pop={value=31.5, method=MLE}
#> add__sd={value=0.7, method=MLE}
#> omega_ka={value=0.774596669241483, method=MLE}
#> omega_cl={value=0.547722557505166, method=MLE}
#> omega_v={value=0.316227766016838, method=MLE}
#> 
#> <MONOLIX>
#> 
#> [TASKS]
#> populationParameters()
#> individualParameters(method = {conditionalMode})
#> fim(method = Linearization)
#> logLikelihood(method = Linearization)
#> plotResult(method = {outputplot, indfits, obspred, residualsscatter, residualsdistribution, parameterdistribution, covariatemodeldiagnosis, randomeffects, covariancemodeldiagnosis, saemresults})
#> 
#> [SETTINGS]
#> GLOBAL:
#> exportpath = 'one.compartment-monolix'
#> POPULATION:
#> exploratoryautostop = no
#> smoothingautostop = no
#> burniniterations = 5
#> exploratoryiterations = 250
#> simulatedannealingiterations = 250
#> smoothingiterations = 200
#> exploratoryalpha = 0
#> exploratoryinterval = 200
#> omegatau = 0.95
#> errormodeltau = 0.95

^{Created on 2024-05-24 with reprex v2.1.0}

[john-harrold]

So 'f$mlxtran' is the mlxtran file. What file extension should 'f$monolixModel" have?

[mattfidler]

It says the file name in the mlxtran file:

file='one.compartment-monolix.txt'

Also the dataset is quite different than rxode2 and NONMEM, you will definitely have to translate to an adm dataset instead of a cmt dataset.

I really suggest using

nlmixr(model, data, "monolix") so it will do all of this for you.

If you want a feature (like not trying to run monolix/nonmem), please let me know as I just started submitting the nlmixr2 stack for PAGE.

[john-harrold]

There are two things in my head :). One is that with ubiquity, I just trranslate models regardless of datasets. So what I want is the ability to take a model I created in ubiquity and say used Matlab to develop (plotting, naive-pooled parameter estimation, etc.) and then say I want to give that to someone in NONMEM format. The underlying mechanics of ubiquity mean I'm already generating the model in rxode2 format. So I've got functions that will go from

ubiqiuty format --> rxode2 output --> NONMEM via babelmixr2.

So there is no dataset here. Just the model. It's about interoperability and trying to give folks something they are familiar with.

The other use-case for this is to assure folks that there is no vendor (I know we're not vendors but you get the idea) lock in for nlmixr2. Sure you can use nlmixr2 but if you decide you want to shift everything to Monolix and only use that, well we can give you models in that format in an automated fashion.

I'm happy using the stub method you provided above, but it might be nice to provide a simple documented way for folks to do that. I don't think it needs to be done before PAGE.

[mattfidler]

I still think nlmixr should be the function. Perhaps without the dataset. This way I use the DRY development principles.

[john-harrold]

If you want to create a better way that's cool. I'm trying not to create work for something that very few people may not use. If you create a method, I'll use that.

[mattfidler]

Based on your use-case you should probably only provide f$monolixModel since you need the data to open a Monolix project.