R- ode函数(deSolve包)：将参数的值更改为时间的函数。

Question

我试图用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.02次单位的步长求解0到48次之间的微分方程。因此，我们指定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

我已经为同样的问题挣扎了一段时间。我试图复制IV输注，然后使用deSolve包，而不是原始模型中使用的RxODE包进行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))})}

我从安德鲁·布克那里得到了Infusion <-线的想法，如https://github.com/andrewhooker/PopED/issues/11所示。我将他的建议从if else修改为ifelse，并在第5次输注结束后加入time <= max(Inf.Times)，否则该模型将继续重新注入。这还允许我在使用deSolve运行模型时使用events表实现额外的非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图相同，为6位数。我能够复制静脉输注(第一组5个峰)，也可以重复PO的剂量(第二组5个峰)。

​

Answer 2

deSolve的大多数求解者使用一个自动的内部时间步长，根据系统的粗糙度或平滑度来调整自己。在模型函数中使用if语句或if()函数并不是一个好主意，原因有两个：(1)时间步长可能不准确；(2)模型函数(即导数)最好是连续可微性，并且避免分步行为，即使在这种情况下求解器是相当稳健的。

deSolve包提供了两种解决问题的方法：“强制函数”和“事件”。两者各有优缺点，但与集成时间步骤相比，如果“事件”(例如注入)的时间非常短，则“事件”尤其有用。

关于这一点的更多信息，可以在deSolve帮助页面?forcings和?events、useR!2017会议的deSolve:强制函数和事件和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添加一个“量”，而“速率”表示在时间间隔中的连续添加。这通常被称为强迫函数。

强迫函数甚至更简单，而且在数值上更好。