我使用 dMod 包在 R 中编写了一个相当复杂的常微分方程药代动力学/药效学模型。我已将其转换为 ODEsensitivity 包可读的方程,因此我可以获得速率常数的 Sobol 灵敏度指数。我可以将速率常数的贡献返回到我在不同隔室中的药物浓度值。
然而,我想要的读数是隔间中药物暴露的持续时间——这可以通过药物超过有效浓度的时间量或曲线下面积来实现。但是,我不明白如何评估任何比药物浓度更复杂的指标的 Sobol 灵敏度。如何返回不同输出的 Sobol 索引?
我写了下面这个简单的玩具版本来模拟两个不同大小的隔间之间的药物转移来演示。ster_C1
下面的代码返回和的每个时间点的敏感度数据框ster_C2
。我宁愿评估一些任意时间点之间的灵敏度ster_C1/ster_C2
或曲线下面积。ster_C1
提前致谢!
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)