R-ode 函数（deSolve 包）：将参数的值更改为时间的函数

Question

我正在尝试使用 deSolve 包中的函数 ode 求解一阶微分方程。问题如下：药物在某些时间（输注时间）以恒定输注速率给药，并以一级速率消除。因此，该过程可以描述为：

if(t %in% Infusion_times){Infusion <- Infusion_rate} else{Infusion <- 0}
    dC <- -Ke*C + Infusion

其中 t 是时间，Infusion_times 是包含给药时间的向量，C 是药物的数量，Ke 是它的消除常数和Infusion是一个变量，当有输液时取输液速率的值，否则取值为0。因此，假设我们想要管理 3 剂，从时间 0、24 和 40 开始，每次输液持续两个小时，假设我们希望 deSolve 每 0.02 个时间单位计算一次答案。例如，我们希望 deSolve 求解 0 到 48 之间时间的微分方程，步长为 0.02 倍。所以我们的向量指定 ode 函数的时间将是：

times <- seq(from = 0, to = 48, by = 0.02)

输液时间由下式给出：

Infusion_times <- c(seq(from = 0, to = 2, by = 0.02), seq(from = 24, to = 26, by = 0.02), 
                    seq(from = 40, to = 42, by = 0.02))

起初，我担心问题可能出在浮点数的比较上。为了防止这种情况，我将两个向量四舍五入到小数点后两位：

times <- round(times, 2)
Infusion_times <- round(times, 2)

所以现在，希望所有 Infusion_times 都包含在 times 向量中：

(sum(Infusion_times %in% times)/length(Infusion_times))*100
[1] 100

如您所见，Infusion_times (100%) 中的所有值都包含在向量 times 中，因此变量 Infusion 应取值 Infusion_rate 在指定的时间。但是，当我们解方程时，它不起作用。让我们证明一下，但首先，让我们指定 ode 函数所需的其他值：

parameters <- c(Ke = 0.5)
amounts <- c(C = 0) #Initial value for drug is 0
Inf_rate <- 5

现在，让我们根据需要编写一个说明变化率的函数：

OneComp <- function(t, amounts, parameters){
  with(as.list(c(amounts, parameters)),{
      if(t %in% Infusion_times){Infuse =Inf_rate} else{Infuse = 0}
      dC <- -Ke*C + Infuse
  list(c(dC))})
}

对于那些不熟悉 deSolve 包的人，函数的参数 t 将说明方程应该积分的次数，amounts 将指定C 的初始值和 parameters 将给出参数的值（在这种情况下，只是 Ke）。现在，让我们求解方程式：

out <- ode(func = OneComp, y = amounts, parms = parameters, times = times)

让我们绘制结果：

library(ggplot2)

ggplot(as.data.frame(out)) + geom_line(aes(x = times, y = C))

如果 Infusion 始终等于 0，这与我们得到的结果完全相同。但是，请注意，如果我们只模拟一次剂量，并且尝试类似的方法，它会起作用：

OneComp <- function(t, amounts, parameters){
      with(as.list(c(amounts, parameters)),{
          if(t < 2){Infuse =Inf_rate} else{Infuse = 0}
          dC <- -Ke*C + Infuse
      list(c(dC))})
    }
out <- ode(func = OneComp, y = amounts, parms = parameters, times = times)
ggplot(as.data.frame(out)) + geom_line(aes(x = times, y = C))

这里，我们让变量Infuse只在时间小于2小时的时候取Inf_rate的值，有效！ 因此，我对这些行为感到非常困惑。这不是更改变量值的问题，也不是浮点数之间比较的问题......知道这些可能是什么吗？谢谢

Answer 1

大多数 deSolve 求解器使用自动内部时间步长，根据系统的粗糙度或平滑度自行调整。在模型函数中使用 if 语句或 if() 函数不是一个好主意，原因有二：(i) 时间步长可能不准确，以及 (2) 模型函数（即导数) 理想情况下应该是 continuous and differentiable 并避免逐步行为，即使求解器在这种情况下非常稳健。

deSolve 包为您的问题提供了两种方法：“强制函数”和“事件”。两者各有利弊，但如果“事件”（例如注射）的时间与积分时间步长相比非常短，则“事件”特别有用。

有关此内容的更多信息，请参见 deSolve 帮助页面 ?forcings 和 ?events，在 deSolve: Forcing functions and Events from the useR!2017 conference, and in the slides 来自 userR!2014。

请检查以下是否适合您：

library("deSolve")

OneComp <- function(t, y, parms){
  with(as.list(c(y, parms)),{
    dC <- -Ke * C
    list(c(dC))
  })
}

eventfunc <- function(t, y, parms) {
  with(as.list(c(y, parms)),{
    C + Inf_rate
  })
}

parms <- c(Ke = 0.5, Inf_rate = 5)

y0 <- c(C = 0)            # Initial value for drug is 0

Infusion_times <- c(seq(from =  0, to =  2, by = 0.02), 
                    seq(from = 24, to = 26, by = 0.02), 
                    seq(from = 40, to = 42, by = 0.02))

# time step can be made bigger to speedup simulation
times <- round(seq(from = 0, to = 48, by = .1) , 2)

# check that all events are in 'times', but no duplicates
# this check is also done by the solver and may print a warning
# to ensure that the user is very careful with this
unique_times <- cleanEventTimes(times, Infusion_times)
times        <- sort(c(unique_times, Infusion_times))

out <- ode(func = OneComp, y = y0, parms = parms, times = times, 
           events = list(func = eventfunc, time = Infusion_times))

plot(out)
rug(Infusion_times)

两条 cleanEventTimes 行是确保模拟命中所有事件时间的一种可能方法。它通常由求解器自动完成，并可能会发出警告，提醒用户对此非常小心。

我使用“基础”图和 rug 来指示注射时间。

我对 Infusion_times 和 Inf_rate 这两个词有些疑惑。在基于事件的方法中，“数量”在离散时间点添加到状态变量 C，而“速率”表示在 时间间隔 内连续添加。这通常称为强制函数。

强制函数会更简单并且在数值上更好。

Answer 2

一段时间以来，我一直在为同样的问题而苦苦挣扎。我试图使用 deSolve 包而不是原始模型中使用的 RxODE 包来复制 IV 输液和 PO 给药。我的解决方案是列出输液时间，然后提取最大值：

tmp <- c()
for (i in seq(from = 1, to = Num.Doses, by = 1)) {
tmp[i] <- (i * Tau) - Tau
tmp2 <- seq(from = 0, to = 2, by = 0.1)
Inf.Times <- round(unlist(lapply(tmp, function(x, Add) x + tmp2)), 3)}

此处，Num.Doses 设置为 5 次静脉输液。 Tau（给药间隔）为 24 小时，0 是输液开始时间，2 是输液结束时间（小时）。

接下来我构建了一个模型，其完整版本是来自 RxODE (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4728294/) 的多室 PKPD 模型：

Model <- function(time, yini, parameters) {
  with(as.list(c(yini, parameters)), {
    Infusion <- ifelse(time %% Tau < Inf.Time & time <= max(Inf.Times), Dose / Inf.Time, 0)
    C1 <- Central / V1
    dDepot <- - (Ka * Depot)
    dCentral <- (Ka * Depot) - (CL * C1) + Infusion
    list(c(dDepot, dCentral))})}

我从 Andrew Booker 那里得到了关于 Infusion <- 行的想法，如 https://github.com/andrewhooker/PopED/issues/11 所示。我将他的建议从 if else 修改为 ifelse，并在第 5 次输注结束后将硬停止与 time <= max(Inf.Times) 合并，否则模型会不断重新输注。这也允许我在使用 deSolve 运行模型时使用事件 table 实施额外的非 IV 给药：

Dose.Events <- data.frame(var = "Depot", time = c(120, 144, 168, 192, 216), value = Dose2, method = "add")
Times <- seq(from = 0, to = Tot.Hours, by = 1)
out <- ode(y = Ini.Con, times = Times, func = Model2, parms = Parms, events = list(data = Dose.Events))

当运行使用完整模型时，输出与使用RxODE中的原始代码几乎相同，更直接和“干净”。根据 AUC 判断，差异最小，sig figs 相同至 6 位。我能够复制 IV 输液（第一组 5 个峰）并重述 PO 剂量（第二组 5 个峰）。

R-ode 函数（deSolve 包）：将参数的值更改为时间的函数

R- ode function (deSolve package): change the value of a parameter as a function of time

r

ode

equation-solving