使用deSolve求解病毒动力学ODE时的浮点误差与阈值处理问题

你当前的代码中，在ODE函数里直接修改i0的局部变量不会影响deSolve求解器维护的状态变量——求解器每次调用你的函数时，都会传入当前的状态值，你修改局部的i0只是临时改变了这个函数内的变量，下一次迭代仍然会用求解器内部的状态值，所以你的条件判断完全起不到作用。

要实现“当数值小于0.001时置为0并继续求解”，需要基于处理后的状态值计算导数，或者直接限制导数的行为，让状态变量能收敛到0。以下是修正后的代码和关键说明：

修正后的代码 #

pars <- c(lambda = 0.272,      # 未感染细胞日产生率
          beta = 0.00027,  # 病毒感染率（每日）
          dt = 0.00136, # 未感染细胞日死亡率
          di = 0.33, # 感染细胞日死亡率
          k = 100, # 游离病毒日产生率
          c = 2, # 病毒自然日死亡率
          beta_treatment = 0,
          drug_start = 50, 
          drug_end = 365)

times <- seq(0, 3650, 1)

startC <- c(t0 = 1000, i0 = 0, v0 = 500) # t: 靶细胞, i: 感染细胞, v: 病毒

ViralDynamics <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    # 先处理感染细胞的状态：低于阈值则置为0
    i_processed <- ifelse(i0 < 0.001, 0, i0)
    
    # 根据是否在治疗期选择感染率参数
    current_beta <- ifelse(t >= drug_start & t <= drug_end, beta_treatment, beta)
    
    # 计算各变量的导数，使用处理后的i_processed
    dtdt <- lambda - current_beta * t0 * v0 - dt * t0
    didt <- current_beta * t0 * v0 - di * i_processed
    # 当i_processed为0时，dvdt直接为 -c*v0，因为没有病毒产生
    dvdt <- k * i_processed - c * v0
    
    return(list(c(dtdt, didt, dvdt)))
  })
}

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

关键修改点 #

• 状态预处理：创建i_processed变量，将低于0.001的i0替换为0，后续所有导数计算都基于这个处理后的值，确保感染和病毒产生过程在i0足够小时停止。
• 代码冗余消除：用ifelse统一处理治疗期和非治疗期的beta参数，避免重复代码块。
• 导数逻辑对齐：当i_processed为0时，dvdt不再有病毒产生项，符合生物学逻辑，同时didt的感染项也会随current_beta调整，治疗期beta_treatment为0时，didt仅保留感染细胞的死亡项（直到i0归零）。

额外优化建议 #

如果希望确保状态变量严格不小于0（避免数值误差导致的负值），可以在函数内添加状态约束：

# 确保状态变量非负（针对数值求解的微小负值）
state <- pmax(state, 0)

不过更推荐通过导数逻辑控制，比如当靶细胞数接近0时，调整dtdt的计算，避免出现负的靶细胞数。

内容的提问来源于stack exchange，提问作者Akshat Arora