How can I calculate the Sobol sensitivities of the area under a curve in R?

96 Views Asked by At

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)
0

There are 0 best solutions below