deSolve ODE 积分错误，我是否使用了错误的函数？

如何解决deSolve ODE 积分错误，我是否使用了错误的函数？

我正在尝试求解一组与生物过程相关的方程。一个方程（约 5 个）用于 C = Co(exp(k1*t)-exp(k2*t) 形式的药代动力学 (PK) 曲线。需要同时求解该方程的导数以及一些酶结合方程和不符合预期的初始结果。故障排除后，意识到如果 k 为负，使用 desolve ode 函数，PK 导数不会自行进行数值积分。我已经尝试了 ode 函数中的所有方法（lsode、lsoda 等），但都没有成功。我试过调整 rtol，它没有解决。

是否有我应该研究的 deSolve ode 函数的替代方法？或者另一种解决这个问题的方法？

下面是带有简化方程的代码来演示问题。 当 k 为负时，积分解与解析结果不匹配。 当 k 为正时，结果符合预期。

第一幅图像，k=0.2 时的结果：当 k 为正时，分析和综合结果匹配

第二个图像，k=-0.2 时的结果：当 k 为负时，积分结果与分析不匹配

library(deSolve)

abi <- function(t,state,parameters) {
    with(as.list(c(state,parameters)),{
        
        dI <- k*exp(k*t)
                list(c(dI))
    })
}

k <- c(-0.2)

times <- seq(0,24,by = 1)

I_analytical <- exp(k*times)

parameters <- c(k)
state <- c(I = 0)

out <- ode(y = state,times = times,func = abi,parms = parameters)

plot(out)
points(I_analytical ~ times)

指出初始条件很容易解决上面的例子，非常有帮助。这是我无法准确积分的方程，我尝试了几个不同的初始条件，但没有真正成功。

library(deSolve)

## Chaos in the atmosphere
CYP <- function(t,{
        #dE <- ksyn - (kdeg * E) + (k2 * EI) - (k1 * E * I)
        #dEI <- (k1 * E * I) - (k2 * EI) + (k4 * EIstar) - (k3 * EI)
        #dEIstar <- (k3 * EI) - (k4 * EIstar)
        #dOcc <- dEI + dEIstar
        dI <- a*tau1*exp(tau1*t) + b*tau2*exp(tau2*t) + c*tau3*exp(tau3*t)
            #list(c(dE,dEI,dEIstar,dOcc,dI))
          list(c(dI))
    })
}

ifit <- c(-0.956144311,0.82619445,0.024520276,-0.913499862,-0.407478829,-0.037174745)
a = ifit[1]
b = ifit[2]
c = ifit[3]
tau1 = ifit[4]
tau2 = ifit[5]
tau3 = ifit[6]


parameters <- c(ksyn = 0.82,kdeg = 0.02,k1 = 2808,k2 = 370.66,k3 = 2.12,k4 = 0.017,a,b,c,tau1,tau2,tau3)

#state <- c(E = 41,EI = 0,EIstar = 0,Occupancy = 0,I = 0.0)
state <- c(I=-0.01)
times <- seq(0,by = .1)
out <- ode(y = state,func = CYP,parms = parameters)

I_analytical <- a*exp(tau1*times) + b*exp(tau2*times) + c*exp(tau3*times)

plot(out)
points(I_analytical ~ times)

Target curve and the ode solution line.

解决方法

初始值应该是

state <- c(I= a + b + c)
#state <- c(I = 1)

,

第一个脚本包含几个问题。最重要的两个是 (1) 模型函数 (abi) 必须包含导数，而不是积分函数，而 (2) 分析积分模型错过了由积分常数产生的 I_0。

假设一阶衰减模型

dI/dt = k I

然后分析积分产生

I_t = I_0 exp(kt)

代码如下：

library(deSolve)

abi <- function(t,state,parameters) {
  with(as.list(c(state,parameters)),{
    # dI <- k*exp(k*t) # original
    dI <-  k * I       # corrected,should be the dervivative
    list(c(dI))
  })
}

k <- -0.2 # simplified,c() was not necessary
times <- seq(0,24,by = 1)

# correction: set I0 to a value > zero
I0 <- 10

# I_analytical <- exp(k*times)    # original
I_analytical <- I0 * exp(k*times) # corrected,multiplied with I0

#state <- c(I = 0) # original
state <- c(I = I0) # corrected
parameters <- c(k = k)

out <- ode(y = state,times = times,func = abi,parms = parameters)

plot(out)
points(I_analytical ~ times)

如果您愿意，可以进一步简化此代码。