R语言deSolve求解ODE时导数数量与初始条件向量长度不匹配问题

报错信息：

Error in checkFunc(Func2, times, y, rho) :  The number of derivatives returned by func() (3006) must equal the length of the initial conditions vector (6)

函数返回的导数大小等于tspan长度加初始条件向量数量。已尝试调整tspan定义位置、修改迭代步长（从0.1改为1），问题仍未解决。

原代码 #

library(deSolve)
tend <- 300
dt <- 1
tspan <- seq(0, tend, by = dt)

# Parameters
rHI <- 0.4 # PI producers
rHN <- 0.4 # X2 free riders
rLI <- rHI * 0.1 # Xr resistant
rLN <- rHI * 0.1
K <- 10e13 # carrying capacity of tumor

delta0 <- 2
rhoS <- 1
rhoI <- 1
lambdaA <- 0.5
lambdaS <- 0.5
beta <- 1
eta <- 1e3
lambdaAS <- 0.02
BCG_dose <- 4e6
t0 <- 0      
tnow <- times[1]


## Initial Conditions
TV <- 1.4 * 10e10
HI0 <- 0.5 * TV
HN0 <- 0.0 * TV
LI0 <- 0.25 * TV
LN0 <- 0.25 * TV
A0 <- 0
S0 <- 0


##Setting Up 
initcond <- c(HI0, HN0, LI0, LN0, A0, S0)
names(initcond) <- c("HI", "HN", "LI", "LN", "A", "S")
##Setting Up BCG dose
interval <- 7
BCG <- rep(0, length(tspan))
BCG[seq(interval, length(tspan), by = interval)] <- BCG_dose
for (i in 2:length(tspan)) {
  if (tspan[i] - tnow >= interval) {
    # update tnow and BCG0 for exponential decay
    BCG0 <- BCG[i-1]
    t0 <- tnow
    tnow <- tspan[i]
    BCG[i] <- BCG0 * exp(-0.3 * (tnow - t0))
  } else {
    BCG[i] <- BCG[i-1]
  }
}
bladder_ODE <- function(t, x, params, BCG) {
  HI <- x[1]
  HN <- x[2]
  LI <- x[3]
  LN <- x[4]
  A <- x[5]
  S <- x[6]
  
  BCG0 <- BCG[which.max(tspan <= t)]
  dBCG <- -BCG0 * exp(-0.3 * (t - t0))
  BCG[which.max(tspan <= t)] <- BCG0 + dBCG
  
  delta <- delta0 * max((A/(S+A) - 0.5), 0)
  
  dHI <- rHI * HI * (1 - (HI + HN + LI + LN) / K) - delta * HI
  
  dHN <- rHN * HN * (1 - (HI + HN + LI + LN) / K)
  
  dLI <- rLI * LI * (1 - (HI + HN + LI + LN) / K)
  
  dLN <- rLN * LN * (1 - (HI + HN + LI + LN) / K)
  
  dA <- rhoS * (HI + LI) + beta * BCG / (eta + BCG) - lambdaA * A
  
  dS <- rhoI * (HN + LN) - lambdaS * S - lambdaAS * A * S
  
  return(list(c(dHI, dHN, dLI, dLN, dA, dS)))
}

out <- ode(y = initcond, times = tspan, func = bladder_ODE)

问题根源 #

1. ODE函数修改全局变量BCG：函数内执行BCG[which.max(tspan <= t)] <- BCG0 + dBCG，直接修改了全局的BCG向量。后续计算dA时使用整个BCG向量而非单个数值，导致dA变成向量，最终返回的导数总长度变为6+3000=3006，触发报错。
2. 参数传递错误：调用ode时未传递BCG参数，且函数定义的params参数未实际使用。
3. 未定义变量报错：原代码中tnow <- times[1]的times未定义，应改为tspan[1]。

修正步骤 #

1. 预先计算BCG浓度：在ODE函数外计算所有时间点的BCG值，避免在函数内修改全局变量。
2. 正确传递参数：将BCG、tspan和模型参数通过ode的参数传递给函数。
3. 确保导数用单个数值计算：在ODE函数中根据当前时间提取对应的BCG单值，而非使用整个向量。

修正后代码 #

library(deSolve)
tend <- 300
dt <- 1
tspan <- seq(0, tend, by = dt)

# 打包所有参数
params <- list(
  rHI = 0.4,
  rHN = 0.4,
  rLI = 0.4 * 0.1,
  rLN = 0.4 * 0.1,
  K = 10e13,
  delta0 = 2,
  rhoS = 1,
  rhoI = 1,
  lambdaA = 0.5,
  lambdaS = 0.5,
  beta = 1,
  eta = 1e3,
  lambdaAS = 0.02,
  BCG_dose = 4e6,
  interval = 7
)

# 初始条件
TV <- 1.4 * 10e10
initcond <- c(
  HI = 0.5 * TV,
  HN = 0.0 * TV,
  LI = 0.25 * TV,
  LN = 0.25 * TV,
  A = 0,
  S = 0
)

# 预先计算所有时间点的BCG浓度
compute_BCG <- function(times, params) {
  BCG <- rep(0, length(times))
  interval <- params$interval
  BCG_dose <- params$BCG_dose
  tnow <- times[1]
  
  # 标记给药时间点
  dose_indices <- seq(interval, length(times), by = interval)
  BCG[dose_indices] <- BCG_dose
  
  for (i in 2:length(times)) {
    if (times[i] - tnow >= interval) {
      tnow <- times[i]
      # 给药后开始衰减
      BCG[i] <- BCG_dose * exp(-0.3 * (times[i] - tnow))
    } else {
      # 非给药时间继续衰减
      BCG[i] <- BCG[i-1] * exp(-0.3 * (times[i] - times[i-1]))
    }
  }
  return(BCG)
}

BCG <- compute_BCG(tspan, params)

# 修正后的ODE函数
bladder_ODE <- function(t, x, params, BCG, tspan) {
  with(as.list(c(x, params)), {
    # 获取当前时间对应的BCG单值
    bcg_val <- BCG[findInterval(t, tspan)]
    
    delta <- delta0 * max((A/(S+A) - 0.5), 0)
    total_tumor <- HI + HN + LI + LN
    
    dHI <- rHI * HI * (1 - total_tumor / K) - delta * HI
    dHN <- rHN * HN * (1 - total_tumor / K)
    dLI <- rLI * LI * (1 - total_tumor / K)
    dLN <- rLN * LN * (1 - total_tumor / K)
    dA <- rhoS * (HI + LI) + beta * bcg_val / (eta + bcg_val) - lambdaA * A
    dS <- rhoI * (HN + LN) - lambdaS * S - lambdaAS * A * S
    
    return(list(c(dHI, dHN, dLI, dLN, dA, dS)))
  })
}

# 调用ode并传递所有参数
out <- ode(y = initcond, times = tspan, func = bladder_ODE, 
           parms = params, BCG = BCG, tspan = tspan)

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