使用 deSolve 和函数时的浮点误差、阈值

Question

我正在尝试构建一个简单的病毒动力学模型，为此我编写了一个对三种不同速率的 ODE 进行建模的函数。由于某种原因，尽管其中一个常微分方程中有一个负项，当我使用

ode

包中的

deSolve

函数求解它时，它实际上从未达到 0，这会弄乱模型，因为它对任何扰动都非常敏感。

当值达到某个值时，我希望函数将其替换为 0，然后继续。我知道我需要在某个地方放置一些 if 语句，但我已经尝试了所有方法，但它不起作用，并且每次运行模拟时模型都保持不变。

这是该功能，以防有帮助：

pars <- c(lambda = 0.272,      # production of uninfected cells per day
          beta = 0.00027,  # rate of infection by virus per day
          dt = 0.00136, # rate of death of uninfected cells per day
          di = 0.33, # rate of death of infected cells per day
          k = 100, # rate of production of free virus per day
          c = 2, # rate of natural virus death rate per day
          beta_treatment = 0,
          drug_start = 50, 
          drug_end = 365)

times <- seq(0, 3650, 1)

startC <- c(t0 = 1000, i0 = 0, v0 = 500) # t = target cell, i = infected cell, v = virus

ViralDynamics <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    if (t >= drug_start & t <= drug_end) {
      dtdt <- lambda - beta_treatment*t0*v0 - dt*t0 # rate of t cell production per day 
      didt <- beta_treatment*t0*v0 - di*i0 # rate of conversion of t cell to infected t cell
      dvdt <- k*i0 - c*v0 # rate of viral production per day 
      if (i0 < 0.001) {
        i0 <- 0
      }
      return(list(c(dtdt, didt, dvdt)))
    } else {
      dtdt <- lambda - beta*t0*v0 - dt*t0 # rate of t cell production per day 
      didt <- beta*t0*v0 - di*i0 # rate of conversion of t cell to infected t cell
      dvdt <- k*i0 - c*v0 # rate of viral production per day 
      if (i0 < 0.001) {
        i0 <- 0
      }
      return(list(c(dtdt, didt, dvdt)))
    }
  })
}

solution <- ode(y = startC, times = times, func = ViralDynamics, parms = pars,
                method = "ode45")

Answer 1

我认为主要原因是在模型函数结束时将

i0

设置为零。 ode 模型函数仅向求解器返回导数，因此这会被忽略，因为求解器负责处理状态。这就是它通常的工作方式，详细信息将是一个单独的主题。

可以通过事件函数或根函数触发的事件函数来将状态变量任意设置为固定值。 deSolve 帮助页面和一个小教程中对此进行了描述。

它需要更多一点的开销，所以我认为对于这里的问题，我们可以使用

if

语句实用地实现它，因为知道微分方程模型中的这种不连续性并不理想。

ViralDynamics <- function(t, state, parameters){
  with(as.list(c(state, parameters)), {
  
    if (t>= drug_start & t<= drug_end) beta <- beta_treatment
    if (i0 < 0.001) i0 <- 0
  
    dtdt <-  lambda - beta * t0 * v0 - dt*t0 
    didt <- beta * t0 * v0 - di*i0
    dvdt <- k*i0 - c*v0 
    return(list(c(dtdt, didt, dvdt), beta=beta))
  })
}

solution <- ode(y = startC, times = times, func = ViralDynamics, parms = pars, method = "lsoda")

plot(solution)

更多建议：

• 与对整个 ODE 系统的两个版本使用

  if

  不同，使用广义公式并仅更改参数

  beta

  会更紧凑且不易出错。

• ode45

  求解器可以工作，但默认的

  lsoda

  更好，即更稳健、更快。
• 可以将额外的“辅助”值添加到返回的列表中以进行调试和绘图，例如

  beta

  。
• 最后，我想知道模拟开始时发生的非常戏剧性的变化。这种行为可能表明不切实际的参数或初始值，但众所周知，病毒速度很快。