使用deSolve模拟生物系统的代码问题排查求助

问题概述 #

基于生物系统假设构建deSolve微分方程模型后，仿真输出与预期趋势不符，以下是针对问题的排查和修正方案：

相关参考图表 #

• 生物系统示意图：
  
• 预期仿真图：
  
• 当前仿真图：
  

当前运行代码 #

#1. 定义初始条件
#2. 定义时间步长
#3. 定义所有状态变量的微分方程
#4. 模拟模型

#----------------------------------------------------------------------------
library(deSolve)

#定义初始条件
states <- c(R=1, 
            Rp=0,
            RS=1,
            RSp=0,
            S=1)

#定义时间步长
times <- c(seq(0, 20, 0.1))

#定义模型的常微分方程
model1 <- function(time, states, parameters) {
  with (as.list(c(states, parameters)), {
    dR = -(R*k1*S) + (Rp*k2) + (Rp*kfeed*RSp)      
    dRp = (R*k1*S) - (Rp*k2) - (Rp*kfeed*RSp)
    dRS = -(RS*k5*Rp) + (RSp*k4)
    dRSp = (RS*k5*Rp) - (RSp*k4) - (Rp*kfeed*RSp)
    dS = 0;
    return(list(c(dR, dRp, dRS, dRSp, dS)))
  })
}

#模拟模型
simModel1 <- function(parameters) {
  return(as.data.frame(ode(y = states, times = times, func = model1, parms = parameters)))
}

pstart1 <- c(k1=1, 
             k2=0.001,
             kfeed=100,
             k4=0.01,
             k5=0.01)

sim1 <- simModel1(pstart1)

plot(sim1$time, sim1$Rp, type = "l", xlab = "time", ylab = "Rp", col = 'violet')

问题排查与修正 #

1. 微分方程反应方向错误 #

对照生物系统示意图，kfeed对应的是RSp分解生成Rp的反馈反应，但当前代码完全搞反了该反应的物质流向：

• 错误地将Rp*kfeed*RSp作为消耗Rp、生成R的项，同时消耗RSp，这与示意图逻辑相悖。
• 正确的反应逻辑：RSp通过kfeed反应转化为Rp，因此：
  • dRSp需减去kfeed*RSp（RSp被消耗）
  • dRp需加上kfeed*RSp（Rp被生成）
  • dR不受该反应影响

修正后的微分方程函数：

model1 <- function(time, states, parameters) {
  with (as.list(c(states, parameters)), {
    dR = -(R*k1*S) + (Rp*k2)      
    dRp = (R*k1*S) - (Rp*k2) + (kfeed*RSp)
    dRS = -(RS*k5*Rp) + (RSp*k4)
    dRSp = (RS*k5*Rp) - (RSp*k4) - (kfeed*RSp)
    dS = 0
    return(list(c(dR, dRp, dRS, dRSp, dS)))
  })
}

2. 参数量级不合理 #

当前参数设置会导致仿真趋势偏离预期：

• kfeed=100过大，RSp会被快速消耗，无法形成有效的反馈
• k2=0.001过小，Rp逆向转化为R的速率几乎可以忽略

建议调整参数为更合理的量级，示例：

pstart1 <- c(k1=1, 
             k2=0.1,
             kfeed=1,
             k4=0.01,
             k5=0.1)

3. 物质守恒验证 #

可以添加代码验证模型的物质守恒性（R+Rp+RS+RSp的总量应保持初始值2），确保方程逻辑自洽：

sim1$total <- sim1$R + sim1$Rp + sim1$RS + sim1$RSp
plot(sim1$time, sim1$total, type="l", xlab="时间", ylab="总物质浓度")

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