I have written a fairly complex ordinary differential equation pharmacokinetic / pharmacodynamic model in R using the dMod package. I’ve converted this into equations readable by the ODEsensitivity package so I can get Sobol sensitivity indices for my rate constants. I can return the contribution of the rate constants to the value of my drug concentration in different compartments.
However, my desired readout is a duration of drug exposure in a compartment – this could be achieved by the amount of time for which the drug is over an effective concentration, or by an area-under-the-curve. However, I don’t understand how I can evaluate a Sobol sensitivity for any metric more complicated than the drug concentration. How can I return the Sobol indices for a different output?
I wrote the following simple toy version to simulate drug transfer between two compartments of different size to demonstrate. The code below returns a data frame of sensitivities at each point in time for ster_C1
and ster_C2
. I would rather evaluate the sensitivity of ster_C1/ster_C2
, or the area under the curve for ster_C1
between some arbitrary timepoints. Thanks in advance!
library(ODEsensitivity)
toymod <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
dchem_B <- -1*(k01*chem_B) +1*(k02*chem_C)*(1/10)
dchem_C <- 1*(k01*chem_B)*(10/1) -1*(k02*chem_C)
return(list( c(dchem_B, dchem_C) ))
})
}
toyparams <- c(
chem_B = 1e-10,
chem_C = 0,
k01 = 200,
k02 = 1000
)
toyparnames <- names(toyparams)
toylower <- toyparams/10
toyupper <- toyparams*10
toyinit <- c(chem_B = 1e-10, chem_C = 0)
toytimes <- seq(5e-05, 0.005, len = 100)
set.seed(59281)
toy_sobol <- ODEsobol(mod = toymod, pars = toyparnames,
state_init = toyinit,
times = toytimes,
n = 500,
rfuncs = "runif",
rargs = paste0("min = ", toylower,
", max = ", toyupper),
sobol_method = "Martinez",
ode_method = "lsoda",
parallel_eval = TRUE,
parallel_eval_ncores = 2)