解决 ODE 时避免负值

Question

我正在尝试对由 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")

Answer 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

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

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