解决 ODE 时避免负值

Posted 2023-02-18

【中文标题】解决 ODE 时避免负值

【英文标题】：avoid negative values when resolving a ODE

【发布时间】：2022-01-11 14:06:46

【问题描述】：

我正在尝试对由 5 个基因组成的网络的行为进行建模，但我遇到的问题是我得到负值，这在生物学上没有意义。

有没有办法将值限制为零？

我在表示图形时设法做到了，但我不知道如何在主方程中使用 ifelse。

非常感谢-1

###################################################
###preliminaries
###################################################

library(deSolve)
library(ggplot2)
library(reshape2)

###################################################
### Initial values
###################################################

values <- c(A = 1,
            B = 1,
            D = 1,
            E = 20,
            R = 1)

###################################################
### Set of constants
###################################################

constants <- c(a = 1.2,
               b = 0.5,
               c = 1.2,
               d = 1.5,
               e = 0.3,
               f = 0.5,
               g = 1.5,
               h = 0.9,
               i = 1.3,
               j = 1.3,
               m = 0.8,
               n = 0.6,
               q = 1,
               t = 0.0075,
               u = 0.0009,
               Pa = 100,
               Pb = 0.05,
               Pd = 0.1,
               Pe = 10)

###################################################
### differential equations
###################################################

Dynamic_Model<-function(t, values, constants) 
  with(as.list(c(values, constants)),
    
    dA <- Pa + a*D - j*A - R
    dB <- Pb + b*A + e*E - m*B 
    dD <- Pd + d*B + f*E - g*A - n*D
    dE <- Pe - h*B + i*E - q*E
    dR <- t*A*B - u*D*E 
    
    list(c(dA, dB, dD, dE, dR))
  )   


###################################################
### time
###################################################

times <- seq(0, 200, by = 0.01)

###################################################
### print ## Ploting
###################################################

out <- ode(y = values, times = times, func = Dynamic_Model, parms = constants)

out2 <- ifelse(out<0, 0, out)

out.df = as.data.frame(out2) 

out.m = melt(out.df, id.vars='time') 
p <- ggplot(out.m, aes(time, value, color = variable)) + geom_point(size=0.5) + ggtitle("Dynamic Model")

【问题讨论】：

下面的帖子展示了如何避免负值：***.com/a/56692927/3677576 并解释了如何做到这一点。

从模型的结构中看不到为什么即使是精确的解决方案也应该留在全正区域。例如，平面A=0 的区域为R>Pa+a*D，其中矢量场指向负值A。对于D 和E 的方程也同样成立，最后一个方程在R=0 平面上有一个非线性分离面。跟踪这一点的最简单方法是使用事件。

是的，除了不知道如何编程之外，我还有一个问题是我真的不懂数学。

【参考方案1】：

我完全同意@Lutz Lehmann 的观点，即负值是模型结构的结果。

方程组允许导数仍然为负，即使状态已经低于零，即状态可以进一步减少。我们没有关于状态是什么的信息，所以下面只是一个技术演示。这里实现了一个无量纲 Monod 型反馈函数fb 作为保护措施。它通常接近一。 km 值应该足够小，以便仅对接近于零的状态值起作用，并且它不应该太小以避免出现数值错误。它可以为每个状态单独制定。其他函数类型也是可能的。

library(deSolve)
library(ggplot2)
library(reshape2)

values <- c(A = 1,
            B = 1,
            D = 1,
            E = 20,
            R = 1)

constants <- c(a = 1.2,
               b = 0.5,
               c = 1.2,
               d = 1.5,
               e = 0.3,
               f = 0.5,
               g = 1.5,
               h = 0.9,
               i = 1.3,
               j = 1.3,
               m = 0.8,
               n = 0.6,
               q = 1,
               t = 0.0075,
               u = 0.0009,
               Pa = 100,
               Pb = 0.05,
               Pd = 0.1,
               Pe = 10,
               km = 0.001)

Dynamic_Model<-function(t, values, constants) 
  with(as.list(c(values, constants)),
    
    fb <- function(x) x / (x+km) # feedback
    
    dA <- (Pa + a*D - j*A - R)         * fb(A)
    dB <- (Pb + b*A + e*E - m*B)       * fb(B)
    dD <- (Pd + d*B + f*E - g*A - n*D) * fb(D)
    dE <- (Pe - h*B + i*E - q*E)       * fb(E)
    dR <- (t*A*B - u*D*E)              * fb(R)
    
    list(c(dA, dB, dD, dE, dR))
  )   


times <- seq(0, 200, by = 0.1)

out <- ode(y = values, times = times, func = Dynamic_Model, parms = constants)
plot(out)

其他提示：

之后删除负值 (out2 <- ifelse(out<0, 0, out)) 是错误的。 移除模型函数中的负值，即

在 main 中使用 ifelse

也可能是错误的，因为它可能导致严重违反质量平衡。

时间步不必非常小。无论如何，它们都会被求解器自动调整。时间步太小会使您的模型变慢，并且您会根据需要获得更多输出。 您的某些参数非常大，因此模型变得非常僵硬。

【讨论】：

谢谢！信息量很大。