deSolve ODE Integration Error，我是不是使用了错误的函数？

Question

我正在尝试求解一组与生物过程相关的方程式。一个方程式（约 5）用于 C = Co(exp(k1*t)-exp(k2*t) 形式的药代动力学 (PK) 曲线。需要同时求解该方程的导数以及一些酶结合方程和不符合预期的初始结果。排除故障后，意识到如果使用 desolve ode 函数 k 为负，则 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

初始值应该是

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

Answer 2

第一个脚本包含几个问题。最重要的两个是 (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)

如果需要，可以进一步简化此代码。

deSolve ODE Integration Error，我是不是使用了错误的函数？

deSolve ODE Integration Error, am I using the wrong function?

r

ode

desolve