knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>")

# Setup the moodels from PopED's
# https://andrewhooker.github.io/PopED/articles/model_def_other_pkgs.html
library(babelmixr2)

library(PopED)

library(mrgsolve)

library(PKPDsim)

library(rxode2)

ff_analytic <- function(model_switch,xt,parameters,poped.db){
  with(as.list(parameters),{
    y=xt
    N = floor(xt/TAU)+1
    f=(DOSE/V)*(KA/(KA - CL/V)) *
      (exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) -
         exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))
    return(list( f=f,poped.db=poped.db))
  })
}


code <- '
$PARAM CL=3.75, V=72.8, KA=0.25
$CMT DEPOT CENT
$ODE
dxdt_DEPOT = -KA*DEPOT;
dxdt_CENT = KA*DEPOT - (CL/V)*CENT;
$TABLE double CP  = CENT/V;
$CAPTURE CP
'

moda <- mrgsolve::mcode("optim", code, atol=1e-8, rtol=1e-8,maxsteps=5000)

ff_ode_mrg <- function(model_switch, xt, p, poped.db){
  times_xt <- drop(xt)
  dose_times <- seq(from=0,to=max(times_xt),by=p[["TAU"]])
  time <- sort(unique(c(0,times_xt,dose_times)))
  is.dose <- time %in% dose_times

  data <-
    tibble::tibble(ID = 1,
                   time = time,
                   amt = ifelse(is.dose,p[["DOSE"]], 0),
                   cmt = ifelse(is.dose, 1, 0),
                   evid = cmt,
                   CL = p[["CL"]], V = p[["V"]], KA = p[["KA"]])

  out <- mrgsolve::mrgsim_q(moda, data=data)

  f <-  out$CP

  f <- f[match(times_xt,out$time)]

  return(list(f=matrix(f,ncol=1),poped.db=poped.db))

}

sfg <- function(x,a,bpop,b,bocc){
  parameters=c(
    KA=bpop[1]*exp(b[1]),
    CL=bpop[2]*exp(b[2]),
    V=bpop[3]*exp(b[3]),
    DOSE=a[1],
    TAU=a[2])
  return( parameters )
}

feps <- function(model_switch,xt,parameters,epsi,poped.db){
  f <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))[[1]]
  y = f*(1+epsi[,1])+epsi[,2]
  return(list(y=y,poped.db=poped.db))
}

poped_db_analytic <- create.poped.database(
  ff_fun =ff_analytic,
  fg_fun =sfg,
  fError_fun=feps,
  bpop=c(KA=0.25,CL=3.75,V=72.8),
  d=c(KA=0.09,CL=0.25^2,V=0.09),
  sigma=c(prop=0.04,add=0.0025),
  m=2,
  groupsize=20,
  xt=c( 1,2,8,240,245),
  minxt=c(0,0,0,240,240),
  maxxt=c(10,10,10,248,248),
  bUseGrouped_xt=1,
  a=cbind(DOSE=c(20,40),TAU=c(24,24)),
  maxa=c(DOSE=200,TAU=24),
  mina=c(DOSE=0,TAU=24))


poped_db_ode_mrg <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_mrg)

pk1cmtoral <- PKPDsim::new_ode_model("pk_1cmt_oral") # take from library
ff_ode_pkpdsim <- function(model_switch, xt, p, poped.db){
    #Set up time points for the ODE
    times_xt <- drop(xt)
    dose_times <- seq(from=0,to=max(times_xt),by=p[["TAU"]])
    times <- sort(unique(c(0,times_xt,dose_times)))
  #
    N = length(dose_times)
    regimen = PKPDsim::new_regimen(amt=p[["DOSE"]],n=N,interval=p[["TAU"]])
    design <- PKPDsim::sim(
      ode = pk1cmtoral,
      parameters = c(CL=p[["CL"]],V=p[["V"]],KA=p[["KA"]]),
      regimen = regimen,
      only_obs = TRUE,
      t_obs = times,
      checks = FALSE,
      return_design = TRUE)
    tmp <- PKPDsim::sim_core(sim_object = design, ode = pk1cmtoral)
    f <- tmp$y
    m_tmp <- match(round(times_xt,digits = 6),tmp[,"t"])
    if(any(is.na(m_tmp))){
      stop("can't find time points in solution\n",
           "try changing the digits argument in the match function")
    }
   #
    f <- f[m_tmp]
    return(list(f = f, poped.db = poped.db))
}

modrx <- rxode2::rxode2({
  d/dt(DEPOT) = -KA*DEPOT;
  d/dt(CENT) = KA*DEPOT - (CL/V)*CENT;
  CP=CENT/V;
})

ff_ode_rx <- function(model_switch, xt, p, poped.db){
  times_xt <- drop(xt)
  et(0,amt=p[["DOSE"]], ii=p[["TAU"]], until=max(times_xt)) %>%
    et(times_xt) -> data

  out <- rxode2::rxSolve(modrx, p, data, atol=1e-8, rtol=1e-8,maxsteps=5000,
                 returnType="data.frame")

  f <-  out$CP[match(times_xt,out$time)]

  return(list(f=matrix(f,ncol=1),poped.db=poped.db))

}

poped_db_ode_rx <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_rx)

poped_db_ode_pkpdsim <- create.poped.database(poped_db_analytic,ff_fun = ff_ode_pkpdsim)

Introduction -- using babelmixr2 with PopED

babelmixr2 now introduces a new method that takes rxode2/nlmixr2 models converts them to a PopED database to help with optimal design.

As in the PopED vignette comparing ODE solvers (and their speeds), this section will:

babelmixr2 ODE solution

The first step of a design using babelmixr2 is to tell babelmixr2 about the design being optimized. This is a bit different than what is done in PopED directly. Below I am using the et() function to create the event table like a typical rxode2 simulation, but it is used to specify the study design:

library(babelmixr2)
library(PopED)

e <- et(amt=1, ii=24, until=250) %>%
  et(time=c(1,2,8,240,245)) %>%
  as.data.frame() %>%
  dplyr::mutate(low=c(NA_real_, 0, 0, 0, 240, 240),
                high=c(NA_real_, 10, 10, 10, 248, 248))


print(e)

PopED/babelmixr2 event table description

Here note that time is the design times for the PopED designs, they can include dosing; only observations are considered the time-points. They become the xt parameter in the PopED database (excluding the doses).

We also build on the structure of the rxode2 event table with simulations. In simulations the sampling windows cause random times to be generated inside the sampling windows. For this reason, the last line of code fixes the times to where we want to have the multiple endpoint design.

Therefore, in this dataset the low becomes minxt and high becomes maxxt.

We chose these because they build on what is already know from nlmixr2 and used in and do not require any extra coding.

Other things you may have to include in your PopED model data frame are:

Getting PopED functions from nlmixr2/rxode2 ui function

Once the design is setup, we need to specify a model. It is easy to specify the model using the nlmixr2/rxode2 function/ui below:

# model
f <- function() {
  ini({
    tKA <- 0.25
    tCL <- 3.75
    tV <- 72.8
    Favail <- fix(0.9)
    eta.ka ~ 0.09
    eta.cl ~ 0.25 ^ 2
    eta.v ~ 0.09
    prop.sd <- sqrt(0.04)
    add.sd <- sqrt(0.0025)
  })
  model({
    ka <- tKA * exp(eta.ka)
    v <- tV * exp(eta.v)
    cl <- tCL * exp(eta.cl)
    d/dt(depot) <- -ka * depot
    d/dt(central) <- ka * depot - cl / v * central
    cp <- central / v
    f(depot) <- DOSE * Favail
    cp ~ add(add.sd) + prop(prop.sd)
  })
}

f <- f() # compile/check nlmixr2/rxode2 model


# Create a PopED database for `nlmixr2`:
poped_db_ode_babelmixr2 <- nlmixr(f, e, "poped",
                                  popedControl(a=list(c(DOSE=20),
                                                      c(DOSE=40)),
                                               maxa=c(DOSE=200),
                                               mina=c(DOSE=0)))

Note when creating a PopED database with a model and a design event table, many of the PopED database components are generated for you.

These are not things that are hidden, but things you can access directly from the model or even from the compiled ui. Much of the other options for optimal design can be specified with the popedControl() function.

This can help you understand what babelmixr2 is doing, we will show what is being added:

PopED's ff_fun from babelmixr2

This is the function that is run to generate the predictions:

# The ff_fun can be retrieved from the ui with f$popedFfFun
f$popedFfFun

Some things to note in this function:

rxode2 models generated from babelmixr2

By describing this, you can also see that there are 2 rxode2 models generated for the PopED database. You can see these inside of the PopED database as well.

The first model uses model times to solve for arbitrary times based on design:

summary(poped_db_ode_babelmixr2$babelmixr2$modelMT)

You can also see the model used for solving scenarios with a number of time points greater than the design specification:

summary(poped_db_ode_babelmixr2$babelmixr2$modelF)

You can also see that the models are identical with the exception of requesting modeled times. You can see the base/core rxode2 model form the UI here:

f$popedRxmodelBase

PopED's fg_fun

babelmixr2 also generates PopEDs fg_fun, which translates covariates and parameters into the parameters required in the ff_fun and used in solving the rxode2 model.

# You can see the PopED fg_fun from the model UI with
# f$popedFgFun:
f$popedFgFun

PopED's error function fError_fun

You can see the babelmixr2 generated error function as well with:

f$popedFErrorFun

One really important note to keep in mind is that PopED works with variances instead of standard deviations (which is a key difference between nlmixr2 and PopED).

This means that the exported model from babelmixr2 operates on variances instead of standard deviations and care must be taken in using these values to not misinterpret the two.

The export also tries to flag this to make it easier to remember.

Other parameters generated by PopED

f$popedBpop # PopED bpop
f$popedNotfixedBpop # PopED notfixed_bpop
f$popedD # PopED d
f$popedNotfixedD # PopED notfixed_d
f$popedCovd # PopED covd
f$popedNotfixedCovd # PopED notfixed_covd
f$popedSigma # PopED sigma (variance is exported, not SD)
f$popedNotfixedSigma # PopED notfixed_sigma

The rest of the parameters are generated in conjunction with the popedControl().

linear compartment models in babelmixr2

You can also specify the models using the linCmt() solutions as below:

f2 <- function() {
  ini({
    tV <- 72.8
    tKA <- 0.25
    tCL <- 3.75
    Favail <- fix(0.9)
    eta.ka ~ 0.09
    eta.cl ~ 0.25 ^ 2
    eta.v ~ 0.09
    prop.sd <- sqrt(0.04)
    add.sd <- fix(sqrt(5e-6))
  })
  model({
    ka <- tKA * exp(eta.ka)
    v <- tV * exp(eta.v)
    cl <- tCL * exp(eta.cl)
    cp <- linCmt()
    f(depot) <- DOSE
    cp ~ add(add.sd) + prop(prop.sd)
  })
}

poped_db_analytic_babelmixr2 <- nlmixr(f, e,
                                       popedControl(a=list(c(DOSE=20),
                                                           c(DOSE=40)),
                                                    maxa=c(DOSE=200),
                                                    mina=c(DOSE=0)))

Comparing method to the speed of other methods

library(ggplot2)
library(microbenchmark)

compare <- microbenchmark(
  evaluate_design(poped_db_analytic),
  evaluate_design(poped_db_analytic_babelmixr2),
  evaluate_design(poped_db_ode_babelmixr2),
  evaluate_design(poped_db_ode_mrg),
  evaluate_design(poped_db_ode_pkpdsim),
  evaluate_design(poped_db_ode_rx),
  times = 100L)


autoplot(compare) + theme_bw()

Note that the babelmixr2 ode solution is the fastest ode solver in this comparison. Among other things, this is because the model is loaded into memory and does not need to be setup each time. (As benchmarks, the mrgsolve, and PKPDsim implementations on the PopED's website are included).

Also to me, the speed of all the tools are reasonable. In my opinion, the benefit of the babelmixr2 interface to PopED is the simplicity of using nlmixr2 / rxode2 functional models or fits directly in PopED without relying on conversions.

The interface is a bit different than the traditional PopED interface, and requires a design data-set as well as a popedControl() to setup a PopED database to run all of the PopED tasks. This is because traditionally nlmixr2 takes a dataset, "estimation" method and controls to change estimation method options.

babelmixr2 adopts the same paradigm of model, data, control to be applied to PopED. This should allow easy translation between the systems. With easier translation, hopefully optimal design in clinical trials will be easier to achieve.

Key notes to keep in mind

Where to find more examples

The examples from PopED have been converted to work with babelmixr2 and are available in the package and on GitHub



nlmixr2/babelmixr documentation built on June 11, 2025, 11:40 a.m.