Help with ODE solver error and parameter dimension issue in PopED model

[Beerahee]

Hi PopED team,
I'm working on a parent-metabolite PK model using PopED and deSolve. I’ve defined a custom ODE system and parameter transformation function, and I’m trying to use plot_model_prediction() to visualize the model.
However, I keep getting the following error:
DLSODA- At T (=R1), too much accuracy requested
for precision of machine.. See TOLSF (=R2)
In above message, R1 = 1e-05, R2 = nan
I’ve tried adjusting rtol and atol in the ode() call, but the issue persists. I also ensured that bpop is a matrix with named columns.
Here’s a link to the full R script I’m using:
Any help or suggestions would be greatly appreciated!
Script
library(PopED)
library(deSolve)

Define the ODE system for parent-metabolite model

PK_ODE <- function(time, state, parameters) {
with(as.list(c(state, parameters)), {
dDepot <- -ka * Depot
dCentral <- ka * Depot - (Cl/V1) * Central - (Q/V1) * Central + (Q/V2) * Peripheral - Kpm * Central
dPeripheral <- (Q/V1) * Central - (Q/V2) * Peripheral
dMetabolite <- Kpm * Central - (Clm/Vm) * Metabolite
list(c(dDepot, dCentral, dPeripheral, dMetabolite))
})
}

Define the model prediction function with relaxed tolerances

ff_fun <- function(model_switch, xt, parameters, poped.db) {
with(as.list(parameters), {
times <- sort(unique(c(0, xt)))
state <- c(Depot = DOSE, Central = 0, Peripheral = 0, Metabolite = 0)
ode_out <- ode(y = state, times = times, func = PK_ODE, parms = parameters, rtol = 1e-4, atol = 1e-6)
ode_df <- as.data.frame(ode_out)
y_parent <- ode_df[match(xt, ode_df$time), "Central"] / V1
y_metab <- ode_df[match(xt, ode_df$time), "Metabolite"] / Vm
y <- cbind(y_parent, y_metab)
return(list(y = y, poped.db = poped.db))
})
}

Define the parameter transformation function

sfg <- function(x, a, bpop, b, bocc) {
parameters <- c(
Tlag = bpop["Tlag"] * exp(b[1]),
ka = bpop["ka"] * exp(b[2]),
V1 = bpop["V1"] * exp(b[3]),
Cl = bpop["Cl"] * exp(b[4]),
Q = bpop["Q"] * exp(b[5]),
V2 = bpop["V2"] * exp(b[6]),
Clm = bpop["Clm"] * exp(b[7]),
Kpm = bpop["Kpm"] * exp(b[8]),
Vm = bpop["Vm"],
DOSE = a[1]
)
return(parameters)
}

Define fixed effects (bpop) as a matrix with named columns

bpop_vals <- matrix(c(
0.355, # Tlag
0.418, # ka
133, # V1
59.8, # Cl
82.9, # Q
958, # V2
139, # Clm
0.664, # Kpm
100 # Vm
), nrow = 1, byrow = TRUE)

colnames(bpop_vals) <- c("Tlag", "ka", "V1", "Cl", "Q", "V2", "Clm", "Kpm", "Vm")

Define random effects (IIV)

d_vals <- diag(c(
0.659^2, # Tlag
0.336^2, # ka
0.498^2, # V1
0.565^2, # Cl
0.629^2, # Q
0.582^2, # V2
0.501^2, # Clm
0.347^2 # Kpm
))

Define residual error (proportional)

sigma_vals <- matrix(c(0.376, 0, 0, 0.388), nrow = 2, byrow = TRUE)

Define sampling times

xt <- c(0.5, 1, 2, 4, 6, 8, 12)

Define dose

a <- 100

Create PopED database

poped.db <- create.poped.database(
ff_fun = ff_fun,
fg_fun = sfg,
fError_fun = feps.add.prop,
bpop = bpop_vals,
d = d_vals,
sigma = sigma_vals,
groupsize = 20,
m = 1,
xt = xt,
minxt = rep(0, length(xt)),
maxxt = rep(24, length(xt)),
model_switch = matrix(1, nrow = 1, ncol = length(xt)),
notfixed_bpop = rep(1, length(bpop_vals)),
a = a,
mina = 0,
maxa = 1000
)

Run model prediction plot

plot_model_prediction(poped.db, model_num_points = 500, facet_scales = "free", PI = TRUE)