Running Monolix

Step 0: What do you need to do to have nlmixr2 run Monolix from a nlmixr2 model

To use Monolix with nlmixr2, you do not need to change your data or your nlmixr2 dataset. babelmixr2 will do the heavy lifting here.

You do need to setup how to run Monolix. If you have setup the lixoftConnectors package from Monolix, no further setup is needed. Instead if you run Monolix from the command line for grid processing (for example) you can figure out the command to run Monolix (it is often useful to use the full command path and set it in the options, ie options("babelmixr2.monolix"="monolix") or use monolixControl(runCommand="monolix"). If needed, I prefer the options() method since you only need to set it once. This could also be a function if you prefer (but I will not cover using the function here).

Step 1: Run a nlmixr2 in Monolix

Lets take the classic warfarin example. The model we use in the nlmixr2 vignettes is:

pk.turnover.emax3 <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)
    ##
    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    prop.err <- 0.1
    pkadd.err <- 0.1
    ##
    temax <- logit(0.8)
    tec50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)
    ##
    eta.emax ~ .5
    eta.ec50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5
    ##
    pdadd.err <- 10
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    emax = expit(temax+eta.emax)
    ec50 =  exp(tec50 + eta.ec50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)
    ##
    DCP = center/v
    PD=1-emax*DCP/(ec50+DCP)
    ##
    effect(0) = e0
    kin = e0*kout
    ##
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center
    d/dt(effect) = kin*PD -kout*effect
    ##
    cp = center / v
    cp ~ prop(prop.err) + add(pkadd.err)
    effect ~ add(pdadd.err) | pca
  })
}

Once monolix is run, you can run the nlmixr2 model using Monolix as if it is new estimation method:

fit <- nlmixr(pk.turnover.emax3, nlmixr2data::warfarin, "monolix",
                              monolixControl(modelName="pk.turnover.emax3"))
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ assuming monolix is running because 'pk.turnover.emax3-monolix.txt' is present
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> ✔ done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 27560
#> ℹ monolix parameter history needs exported charts, please export charts

This fit issues an informational tidbit -

• monolix parameter history needs exported charts, please export charts

This will automatically be generated as well when lixoftConnectors package is generated and you have a recent version of Monolix. If you don’t have that information then the important parameter history plots will not be imported and you cannot see those plots.

Just like with the NONMEM translation, the monolixControl() has modelName which helps control the output directory of Monolix (if not specified babelmixr2 tries to guess based on the model name based on the input).

Printing this out this nlmixr2 fit you see:

fit

$$\begin{align*} {ktr} & = \exp\left({tktr}+{eta.ktr}\right) \\ {ka} & = \exp\left({tka}+{eta.ka}\right) \\ {cl} & = \exp\left({tcl}+{eta.cl}\right) \\ {v} & = \exp\left({tv}+{eta.v}\right) \\ {emax} & = expit({temax}+{eta.emax}, {0}, {1}) \\ {ec50} & = \exp\left({tec50}+{eta.ec50}\right) \\ {kout} & = \exp\left({tkout}+{eta.kout}\right) \\ {e0} & = \exp\left({te0}+{eta.e0}\right) \\ {DCP} & = \frac{{center}}{{v}} \\ {PD} & = {1}-\frac{{emax} {\times} {DCP}}{\left({ec50}+{DCP}\right)} \\ effect({0}) & = {e0} \\ {kin} & = {e0} {\times} {kout} \\ \frac{d \: depot}{dt} & = -{ktr} {\times} {depot} \\ \frac{d \: gut}{dt} & = {ktr} {\times} {depot}-{ka} {\times} {gut} \\ \frac{d \: center}{dt} & = {ka} {\times} {gut}-\frac{{cl}}{{v}} {\times} {center} \\ \frac{d \: effect}{dt} & = {kin} {\times} {PD}-{kout} {\times} {effect} \\ {cp} & = \frac{{center}}{{v}} \\ {cp} & \sim prop({prop.err})+add({pkadd.err}) \\ {effect} & \sim add({pdadd.err}){\lor}{pca} \end{align*}$$

Of particular interest is the comparison between Monolix predictions and nlmixr predictions. In this case, I believe that these also imply the models are predicting the same thing. Note that the model predictions are not as close as they were with NONMEM because Monolix does not use the lsoda ODE solver. Hence this small deviation is expected, but still gives a validated Monolix model.

Optional Step 2: Add conditional weighted residuals/focei objf to Monolix

As in the case of NONMEM, this gives some things that are not available to Monolix, like adding conditional weighted residuals:

fit <- addCwres(fit)
#> → Calculating residuals/tables
#> ✔ done

Which will add nlmixr’s CWRES as well as adding the nlmixr2 FOCEi objective function

Because you now have an objective function compared based on the same assumptions, you could compare the performance of Monolix and NONMEM based on objective function.

To be fair, objective function values must always be used with caution. How the model performs and predicts the data is far more valuable.

Optional Step 3: Use nlmixr2 for vpc, reporting, etc.

Also since it is a nlmixr2 object it would be easy to perform a VPC too:

v1s <- vpcPlot(fit, show=list(obs_dv=TRUE), scales="free_y") +
  ylab("Warfarin Cp [mg/L] or PCA") +
  xlab("Time [h]")
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> ℹ The deprecated feature was likely used in the vpc package.
#>   Please report the issue to the authors.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

v2s <- vpcPlot(fit, show=list(obs_dv=TRUE), pred_corr = TRUE, scales="free_y") +
  ylab("Prediction Corrected Warfarin Cp [mg/L] or PCA") +
  xlab("Time [h]")

v1s

v2s

Notes about Monolix data translation

The input dataset expected to be compatible with rxode2 or nlmixr2. This dataset is then converted to Monolix format:

• The combination of CMT and Dose type creates a unique ADM variable.

• The ADM definition is saved in the monolix model file

• babelmixr2 creates a macro describing the compartment, ie compartment(cmt=#, amount=stateName)

• babelmixr2 also creates a macro for each type of dosing:

  • Bolus/infusion uses depot() and adds modeled lag time (Tlag) or bioavailability (p) if specified

  • Modeled rate uses depot() with Tk0=amtDose/rate. babelmixr2 also adds modeled lag time (Tlag) or bioavailability (p) if specified

  • Modeled duration uses depot() with Tk0=dur, also add adds modeled lag time (Tlag) or bioavailability (p) if specified Turning off a compartment uses empty macro