求解 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)) 是错误的。
去除模型函数中的负值，即

use the ifelse in the main

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

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

求解 ODE 时避免负值

avoid negative values when resolving a ODE

r

if-statement

ode

differential-equations