What is the best way to set the initial conditions for different individuals.

[john-harrold]

I have a situation where I need to start a simulation at a different initial condition for each subject. In the example below I'm doing that with the events table (arbitrarily I'm setting the amount in the central compartment to a different value). I'm curious if this is bad way to do this :). Is there something I'm not really thinking about here?

rm(list=ls())
library(nlmixr2lib)
library(rxode2)


mymod = function() {
  description <- "Two compartment PK model with linear clearance using differential equations"
  ini({
    lka <- 0.45 ; label("Absorption rate (Ka)")
    lcl <- 1 ; label("Clearance (CL)")
    lvc <- 3 ; label("Central volume of distribution (V)")
    lvp <- 5 ; label("Peripheral volume of distribution (Vp)")
    lq  <- 0.1 ; label("Intercompartmental clearance (Q)")
    propSd <- 0.5 ; label("Proportional residual error (fraction)")
  })
  model({
    ka <- exp(lka)
    cl <- exp(lcl)
    vc <- exp(lvc)
    vp <- exp(lvp)
    q  <- exp(lq)

    kel <- cl/vc
    k12 <- q/vc
    k21 <- q/vp

    d/dt(depot) <- -ka*depot
    d/dt(central) <-  ka*depot - kel*central - k12*central + k21*peripheral1
    d/dt(peripheral1) <- k12*central - k21*peripheral1
    Cc <- central / vc

    Cc ~ prop(propSd)
  })
}

mymod = addEta(mymod, c("lka", 'lcl', 'lq', 'lvp', 'lvc'))

amts = c(20, 50, 75)

ev = NULL
ot = seq(0,50, .1)
for(idx  in 1:length(amts)){
  evtmp = et(ot, id=idx)                             |>
          et(amt=amts[idx], time=0, cmt=central, evid=4, id=idx) |>
          et(amt=100, time=2, cmt=central, id=idx)

  if(is.null(ev)){
    ev = evtmp
  } else {
    ev = etRbind(ev, evtmp)
  }
}

sim = rxSolve(mymod, events = ev)

plot(sim, central)

[mattfidler]

I don't think there is anything wrong with this approach, though you could simulate different initial conditions in this case by changing bioavailability f(depot) which would could use the same event table. The same can be said about depot(0).

One note, lapply is a bit faster to create such event tables:

rm(list=ls())
library(nlmixr2lib)
library(rxode2)
#> rxode2 2.1.1 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`


mymod = function() {
  description <- "Two compartment PK model with linear clearance using differential equations"
  ini({
    lka <- 0.45 ; label("Absorption rate (Ka)")
    lcl <- 1 ; label("Clearance (CL)")
    lvc <- 3 ; label("Central volume of distribution (V)")
    lvp <- 5 ; label("Peripheral volume of distribution (Vp)")
    lq  <- 0.1 ; label("Intercompartmental clearance (Q)")
    propSd <- 0.5 ; label("Proportional residual error (fraction)")
  })
  model({
    ka <- exp(lka)
    cl <- exp(lcl)
    vc <- exp(lvc)
    vp <- exp(lvp)
    q  <- exp(lq)

    kel <- cl/vc
    k12 <- q/vc
    k21 <- q/vp

    d/dt(depot) <- -ka*depot
    d/dt(central) <-  ka*depot - kel*central - k12*central + k21*peripheral1
    d/dt(peripheral1) <- k12*central - k21*peripheral1
    Cc <- central / vc

    Cc ~ prop(propSd)
  })
}

mymod = addEta(mymod, c("lka", 'lcl', 'lq', 'lvp', 'lvc'))
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ promote `etalka` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalka` to `0.1`
#> ℹ promote `etalcl` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalcl` to `0.1`
#> ℹ promote `etalq` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalq` to `0.1`
#> ℹ promote `etalvp` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalvp` to `0.1`
#> ℹ promote `etalvc` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalvc` to `0.1`

amts = c(20, 50, 75)

ev = NULL
ot = seq(0,50, .1)
f1 <- function() {
  do.call(etRbind, lapply(1:length(amts), function(idx) {
    evtmp = et(ot, id=idx)                             |>
      et(amt=amts[idx], time=0, cmt=central, evid=4, id=idx) |>
      et(amt=100, time=2, cmt=central, id=idx)
  }))
}

f2 <- function() {
  for(idx  in 1:length(amts)){
    evtmp = et(ot, id=idx)                             |>
      et(amt=amts[idx], time=0, cmt=central, evid=4, id=idx) |>
      et(amt=100, time=2, cmt=central, id=idx)

    if(is.null(ev)){
      ev = evtmp
    } else {
      ev = etRbind(ev, evtmp)
    }
  }
  ev
}

sim = rxSolve(mymod, events = f1())
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’

plot(sim, central)

^{Created on 2023-12-18 with reprex v2.0.2}