Working with rxode2, nlmixr2, NONMEM, and
Monolix
Say you have a model in ubiquity and you want to use it with nlmixr2.
Or perhaps you want to try out the simulation engine of
rxode2, or you want to hand the model off to someone using
Monolix. This should provide you with a way to easily
convert your model into different formats in an automated fashion.
Getting your model in rxode2 format
As an example consider the system file for the two compartment model
(system_2cmt.txt below). When you build the system an
rxode2 output target will be created. You can use the
system_fetch_template() function to access this file. Below
I’m placing the script into the temporary directory. This is a script
that will generate an rxode2 model function and store the
function in the object my_model. This is the default input
for nlmixr2. You can then use this function and the piping
methodology to do things like add IIV terms to parameters. So you can
use ubiquity to build your structural model and then
piping, rxode2, and nlmixr2 for analysis.
cfg = build_system(system_file = sf_2cmt_full,
temporary_directory = file.path(tempdir(), "transient"))
fr = system_fetch_template(cfg,
template = "nlmixr2",
output_directory = tempdir(),
overwrite = TRUE)
library(rxode2)
source(file.path(tempdir(), "system_nlmixr2.R"))You can see the contents of system_nlmixr2.R and what
the my_model object looks like below. Note that to include
the output with ubiquity you will need to have some error
model. So along with the <O> Descriptor you will need
to define some variance parameters (<VP>) and also an
output error model (<OE:?>).
Getting your model in NONMEM format
Once you have your model in rxode2 format you can
convert it further into NONMEM format with the babelmixr2
package (via the ruminate package). To do this you need at
least one inter-individual variability term. You can do this in the
system file with the <IIV:?> delimiter. In this
example we will do that by using model piping. Here we specify that the
typical value of Vc, the parameter TV_Vc,
should have an error that is lognormally distributed. This will add the
IIV term eta_Vc.
Next we just need to load the ruminate package
(>=0.2.4) and use the rx2other() function to create a
NONMEM control stream:
Getting your model in Monolix format
You can also get the model in Monolix by specifying
monolix as the output type:
my_model_mlx = rx2other(my_model, out_type="monolix")Models in different formats
Contents of system_2cmt.txt file:
# Two compartment model with absorption compartment (At). The central
# compartment (Cc) had a volume (Vc), and the tissue compartment (Cp) has a
# volume (Vt). The system is parameterized in terms of macro constants (Q and
# CL)
#
# While the the default dosing is a 1 mg IV dose into the central compartment,
# the model is written to accept dosing into an absorption compartment and
# through continuous IV infusion.
# _________
# | |
# | At |
# | |
# |_________|
# |
# | ka, fb
# |
# V
# _________ _________
# | | Q | |
# | Cc |------>| Cp |
# | Vc |<------| Vt |
# |_________| |_________|
# |
# | CL
# |
# V
#
# System Units:
#
# mass [=] mg
# volume [=] ml
# concentration [=] mg/ml
# time [=] hr
# #-------------#
# | Parameters |
# #-------------#
#
# System parameters
# name value lower upper units editable grouping
# bound bound
<P> Vc 1.0 eps Inf ml yes System
<P> Vt 1.0 eps Inf ml yes System
<P> CL 1.0 eps Inf ml/hr yes System
<P> Q 1.0 eps Inf ml/hr yes System
<P> ka 1.0 eps Inf 1/hr yes System
<P> fb 1.0 eps Inf -- yes System
# Bolus Events
# ------------
# times/events state values scale units
<B:times>; [ 0 ]; 1; hours
<B:events>; Cc; [1.0 ]; 1/Vc; mg
<B:events>; At; [0.0 ]; 1; mg
# Infusion Rates
# ------------
# name time/levels values scale units
<R:Dinf>; times; [0]; 1; hours
<R:Dinf>; levels; [0]; 1; mg/hour
<ODE:At> -ka*At
<ODE:Cc> ka*At*fb/Vc - CL/Vc*Cc -Q*(Cc-Cp)/Vc + Dinf/Vc
<ODE:Cp> +Q*(Cc-Cp)/Vt
<O> Cc_mg_ml = Cc
<VP> prop_err 0.1 eps inf -- yes Variance
<VP> add_err 0.1 eps inf ng/ml yes Variance
<OE:Cc_mg_ml> add=add_err; prop=prop_err
<TS:hours> 1.0
<TS:days> 1.0/24Contents of system_nlmixr2.R file produced by
system_fetch_template():
#
#
#
my_model <- function() {
ini({
# Typical Value of System Parameters
TV_Vc = c(.Machine$double.eps, 1.0, Inf)
TV_Vt = c(.Machine$double.eps, 1.0, Inf)
TV_CL = c(.Machine$double.eps, 1.0, Inf)
TV_Q = c(.Machine$double.eps, 1.0, Inf)
TV_ka = c(.Machine$double.eps, 1.0, Inf)
TV_fb = c(.Machine$double.eps, 1.0, Inf)
# Between-subject variability:
# Error model parameters
prop_err = c(.Machine$double.eps, 0.1, .Machine$double.xmax)
add_err = c(.Machine$double.eps, 0.1, .Machine$double.xmax)
})
model({
# System Parameters
Vc = TV_Vc
Vt = TV_Vt
CL = TV_CL
Q = TV_Q
ka = TV_ka
fb = TV_fb
# Placeholders for infusion rates.
# These should be defined in the dataset.
# time scale: 1 units: (hours)
# mass scale: 1 units: (mg/hour)
Dinf = 0.0
# Defining ODEs
d/dt(At) = (-ka*At)
d/dt(Cc) = (ka*At*fb/Vc - CL/Vc*Cc -Q*(Cc-Cp)/Vc + Dinf/Vc)
d/dt(Cp) = (+Q*(Cc-Cp)/Vt)
# Outputs and error models
Cc_mg_ml = Cc
Cc_mg_ml ~ add(add_err) + prop(prop_err)
})
}The my_model object:
[36m──
[39m
[1mrxode2-based free-form 3-cmt ODE model
[22m
[36m──────────────────────────────────────
[39m
──
[1m
[1mInitalization:
[1m
[22m ──
[1mFixed Effects
[22m (
[1m
[34m$theta
[39m
[22m):
TV_Vc TV_Vt TV_CL TV_Q TV_ka TV_fb prop_err add_err
1.0 1.0 1.0 1.0 1.0 1.0 0.1 0.1
[1mOmega
[22m (
[1m
[34m$omega
[39m
[22m):
eta_Vc
eta_Vc 1
[1m
States
[22m (
[1m
[34m$state
[39m
[22m or
[1m
[34m$stateDf
[39m
[22m):
Compartment Number Compartment Name
1 1 At
2 2 Cc
3 3 Cp
──
[1m
[1m
[1m
[1mμ-referencing
[1m
[1m (
[1m
[1m
[34m$muRefTable
[39m
[1m
[1m):
[1m
[22m ──
theta eta level
1 TV_Vc eta_Vc id
──
[1m
[1mModel (Normalized Syntax):
[1m
[22m ──
function() {
ini({
TV_Vc <- c(2.22044604925031e-16, 1)
TV_Vt <- c(2.22044604925031e-16, 1)
TV_CL <- c(2.22044604925031e-16, 1)
TV_Q <- c(2.22044604925031e-16, 1)
TV_ka <- c(2.22044604925031e-16, 1)
TV_fb <- c(2.22044604925031e-16, 1)
prop_err <- c(2.22044604925031e-16, 0.1)
add_err <- c(2.22044604925031e-16, 0.1)
eta_Vc ~ 1
})
model({
Vc <- exp(TV_Vc + eta_Vc)
Vt = TV_Vt
CL = TV_CL
Q = TV_Q
ka = TV_ka
fb = TV_fb
Dinf = 0
d/dt(At) = (-ka * At)
d/dt(Cc) = (ka * At * fb/Vc - CL/Vc * Cc - Q * (Cc -
Cp)/Vc + Dinf/Vc)
d/dt(Cp) = (+Q * (Cc - Cp)/Vt)
Cc_mg_ml = Cc
Cc_mg_ml ~ add(add_err) + prop(prop_err)
})
}The file containing the control stream is found here:
my_model_nm$files$ctl$fn_full.
$PROBLEM translated from babelmixr2
; comments show mu referenced model in ui$getSplitMuModel
$DATA my_model.csv IGNORE=@
$INPUT ID TIME EVID AMT DV CMT RXROW
$SUBROUTINES ADVAN13 TOL=6 ATOL=12 SSTOL=6 SSATOL=12
$MODEL NCOMPARTMENTS=3
COMP(AT, DEFDOSE) ; At
COMP(CC) ; Cc
COMP(CP) ; Cp
$PK
MU_1=THETA(1)
VC=DEXP(MU_1+ETA(1)) ; Vc <- exp(TV_Vc)
VT=THETA(2) ; Vt <- TV_Vt
CL=THETA(3) ; CL <- TV_CL
RXR1=THETA(4) ; Q <- TV_Q
KA=THETA(5) ; ka <- TV_ka
FB=THETA(6) ; fb <- TV_fb
$DES
DINF=0 ; Dinf = 0
DADT(1) = (- KA*A(1)) ; d/dt(At) = (-ka * At)
DADT(2) = (KA*A(1)*FB/VC-CL/VC*A(2)-RXR1*(A(2)-A(3))/VC+DINF/VC) ; d/dt(Cc) = (ka * At * fb/Vc - CL/Vc * Cc - Q * (Cc - Cp)/Vc + Dinf/Vc)
DADT(3) = (RXR1*(A(2)-A(3))/VT) ; d/dt(Cp) = (Q * (Cc - Cp)/Vt)
CC_MG_ML=A(2) ; Cc_mg_ml = Cc
$ERROR
;Redefine LHS in $DES by prefixing with on RXE_ for $ERROR
RXE_DINF=0 ; Dinf = 0
RXE_CC_MG_ML=A(2) ; Cc_mg_ml = Cc
RX_PF1=RXE_CC_MG_ML ; rx_pf1 ~ Cc_mg_ml
; Write out expressions for ipred and w
RX_IP1 = RX_PF1
RX_P1 = RX_IP1
W1=DSQRT((THETA(8))**2+(RX_PF1)**2*(THETA(7))**2) ; W1 ~ sqrt((add_err)^2 + (rx_pred_f_)^2 * (prop_err)^2)
IF (W1 .EQ. 0.0) W1 = 1
IPRED = RX_IP1
W = W1
Y = IPRED + W*EPS(1)
$THETA (2.2204e-16, 1 ) ; 1 - TV_Vc
(2.2204e-16, 1 ) ; 2 - TV_Vt
(2.2204e-16, 1 ) ; 3 - TV_CL
(2.2204e-16, 1 ) ; 4 - TV_Q
(2.2204e-16, 1 ) ; 5 - TV_ka
(2.2204e-16, 1 ) ; 6 - TV_fb
(2.2204e-16, 0.1) ; 7 - prop_err
(2.2204e-16, 0.1) ; 8 - add_err
$OMEGA 1 ; eta_Vc
$SIGMA 1 FIX
$ESTIMATION METHOD=1 INTER MAXEVALS=100000 SIGDIG=3 SIGL=12 PRINT=1 NOABORT
$COVARIANCE
$TABLE ID ETAS(1:LAST) OBJI FIRSTONLY ONEHEADER NOPRINT
FORMAT=s1PE17.9 NOAPPEND FILE=my_model.eta
$TABLE ID TIME IPRED PRED RXROW ONEHEADER NOPRINT
FORMAT=s1PE17.9 NOAPPEND FILE=my_model.predThe file containing the mlxtran output is found here:
my_model_mlx$files$mlxtran$fn_full
<DATAFILE>
[FILEINFO]
file='my_model-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_Cc_mg_ml, type=continuous}
ADM = {use=administration}
<MODEL>
[INDIVIDUAL]
input={Vc_pop, omega_Vc, Vt_pop, CL_pop, Q_pop, ka_pop, fb_pop}
DEFINITION:
Vc = {distribution=logNormal, typical=Vc_pop, sd=omega_Vc}
Vt = {distribution=normal, typical=Vt_pop, no-variability}
CL = {distribution=normal, typical=CL_pop, no-variability}
Q = {distribution=normal, typical=Q_pop, no-variability}
ka = {distribution=normal, typical=ka_pop, no-variability}
fb = {distribution=normal, typical=fb_pop, no-variability}
[LONGITUDINAL]
input={add_err, prop_err}
file='my_model-monolix.txt'
DEFINITION:
rx_prd_Cc_mg_ml={distribution = normal, prediction = rx_pred_Cc_mg_ml, errorModel=combined2(add_err,prop_err)}
<FIT>
data={rx_prd_Cc_mg_ml}
model={rx_prd_Cc_mg_ml}
<PARAMETER>
Vc_pop={value=2.71828182845905, method=MLE}
Vt_pop={value=1, method=MLE}
CL_pop={value=1, method=MLE}
Q_pop={value=1, method=MLE}
ka_pop={value=1, method=MLE}
fb_pop={value=1, method=MLE}
prop_err={value=0.1, method=MLE}
add_err={value=0.1, method=MLE}
omega_Vc={value=1, method=MLE}
<MONOLIX>
[TASKS]
populationParameters()
individualParameters(method = {conditionalMode})
fim(method = StochasticApproximation)
logLikelihood(method = ImportanceSampling)
plotResult(method = {outputplot, indfits, obspred, residualsscatter, residualsdistribution, parameterdistribution, covariatemodeldiagnosis, randomeffects, covariancemodeldiagnosis, saemresults})
[SETTINGS]
GLOBAL:
exportpath = 'my_model-monolix'
POPULATION:
exploratoryautostop = no
smoothingautostop = no
burniniterations = 5
exploratoryiterations = 250
simulatedannealingiterations = 250
smoothingiterations = 200
exploratoryalpha = 0
exploratoryinterval = 200
omegatau = 0.95
errormodeltau = 0.95