使用R语言deSolve包dede求解时滞SEIRVB ODE的问题排查

核心问题：未向dede告知延迟参数 #

deSolve的dede函数需要明确知晓模型中存在的时间延迟长度，否则求解器会直接忽略滞后项，按普通常微分方程求解，导致修改tau后结果无变化。此外，全局变量的使用也可能引发潜在的作用域问题。

具体修复步骤 #

1. 将所有参数（包括tau）放入parms列表 #

避免依赖全局变量，把模型用到的所有参数打包进parms列表，确保函数能正确获取参数：

# 参数定义改为列表形式
parms <- list(
  Lambda = 0.5,      
  beta1 = 1.13921549,    
  beta2 = 2.68228986,      
  beta3 = 2.9627036,      
  beta4 = 1.19686956,     
  beta5 = 0.5,    
  beta6 = 1.34496108,     
  beta7 = 3.40332936,   
  beta8 = 1.48249240,    
  beta9 = 0.04681161,     
  beta10 = 2.32555864,    
  alpha = 0.5,    
  kappa = 0.02933240,    
  d = 0.0047876,          
  r = 0.09871,        
  tau = 7  # 时间延迟
)

2. 修改模型函数，从parms中读取参数 #

调整模型函数，使用with(parms, {})包裹逻辑，不再依赖全局变量：

seirvb_model <- function(t, y, parms) {
  with(parms, {
    # 获取滞后的I值
    if (t < tau) {
      Ilag <- initial_conditions["I"]
    } else {
      # lagvalue第二个参数为变量索引，y[3]对应I，索引正确
      Ilag <- lagvalue(t - tau, 3)
    }
    
    # 计算各状态变量的导数
    dS <- Lambda - beta1 * y[1] * y[2] - (beta3 + alpha) * y[1]
    dE <- beta4 * y[5] * y[2] + beta6 * y[6] * y[2] + beta1 * y[1] * y[2] - beta2 * y[2] * y[3] - (kappa + alpha) * y[2]
    dI <- beta5 * y[5] * Ilag + beta7 * y[6] * y[3] + beta2 * y[2] * y[3] - (d + alpha + r) * y[3]
    dR <- beta9 * y[5] + beta8 * y[6] + r * y[3] + kappa * y[2] - alpha * y[4]
    dV <- beta3 * y[1] - (beta9 + alpha + beta10) * y[5] - beta4 * y[5] * y[2] - beta5 * y[5] * Ilag
    dB <- beta10 * y[5] - beta6 * y[6] * y[2] - beta7 * y[6] * y[3] - (beta8 + alpha) * y[6]
    
    return(list(c(dS, dE, dI, dR, dV, dB)))
  })
}

3. 调用dede时传递delays参数 #

必须通过delays参数告知求解器模型中的延迟时间，这样求解器才会处理滞后项：

# 求解模型，添加delays参数
ode_output <- dede(
  y = initial_conditions, 
  times = times, 
  func = seirvb_model, 
  parms = parms,
  delays = parms$tau  # 传递延迟参数
)

4. 验证修改效果 #

修改parms$tau的值（比如改为14），重新运行代码，观察输出结果的变化，此时时滞应该能正常发挥作用。

完整修复后的代码 #

library(deSolve)

# 参数定义（打包为列表）
parms <- list(
  Lambda = 0.5,      
  beta1 = 1.13921549,    
  beta2 = 2.68228986,      
  beta3 = 2.9627036,      
  beta4 = 1.19686956,     
  beta5 = 0.5,    
  beta6 = 1.34496108,     
  beta7 = 3.40332936,   
  beta8 = 1.48249240,    
  beta9 = 0.04681161,     
  beta10 = 2.32555864,    
  alpha = 0.5,    
  kappa = 0.02933240,    
  d = 0.0047876,          
  r = 0.09871,        
  tau = 7  # 时间延迟
)

# 带时滞的SEIRVB模型函数
seirvb_model <- function(t, y, parms) {
  with(parms, {
    # 获取滞后的I值
    if (t < tau) {
      Ilag <- initial_conditions["I"]
    } else {
      Ilag <- lagvalue(t - tau, 3)  # y[3]对应I变量
    }
    
    # 计算各状态变量导数
    dS <- Lambda - beta1 * y[1] * y[2] - (beta3 + alpha) * y[1]
    dE <- beta4 * y[5] * y[2] + beta6 * y[6] * y[2] + beta1 * y[1] * y[2] - beta2 * y[2] * y[3] - (kappa + alpha) * y[2]
    dI <- beta5 * y[5] * Ilag + beta7 * y[6] * y[3] + beta2 * y[2] * y[3] - (d + alpha + r) * y[3]
    dR <- beta9 * y[5] + beta8 * y[6] + r * y[3] + kappa * y[2] - alpha * y[4]
    dV <- beta3 * y[1] - (beta9 + alpha + beta10) * y[5] - beta4 * y[5] * y[2] - beta5 * y[5] * Ilag
    dB <- beta10 * y[5] - beta6 * y[6] * y[2] - beta7 * y[6] * y[3] - (beta8 + alpha) * y[6]
    
    return(list(c(dS, dE, dI, dR, dV, dB)))
  })
}

# 初始条件
initial_conditions <- c(
  S=0.3,
  E=0.1,
  I=0.006,
  R=0,
  V=0,
  B=0
)

# 时间向量
times <- seq(0, 120, by = 1)

# 求解模型（传递delays参数）
ode_output <- dede(
  y = initial_conditions, 
  times = times, 
  func = seirvb_model, 
  parms = parms,
  delays = parms$tau
)

# 绘图
plot(ode_output, xlab = "Time", ylab = "Population", main = "SEIRVB Model for COVID-19")

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