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

Question

我正在尝试求解一组与生物过程相关的方程。一个方程（约 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, state, parameters) {
    with(as.list(c(state, parameters)), {
        #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, 24, by = .1)
out <- ode(y = state, times = times, 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.

Answer 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)

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

Answer 2

初始值应该是

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