用零替换模型(ODE系统)中的负值 [英] Replacing negative values in a model (system of ODEs) with zero

问题描述

我目前正在使用deSolve求解常微分方程组，并且想知道是否有任何方法可以防止微分变量值降至零以下.我还看过其他一些有关在向量，数据帧等中将负值设置为零的文章，但是由于这是一个生物学模型(T细胞计数变为负数没有意义)，需要从一开始就阻止它发生，以便这些值不会使结果产生偏差，而不仅仅是替换最终输出中的负数.

I'm currently working on solving a system of ordinary differential equations using deSolve, and was wondering if there's any way of preventing differential variable values from going below zero. I've seen a few other posts about setting negative values to zero in a vector, data frame, etc., but since this is a biological model (and it doesn't make sense for a T cell count to go negative), I need to stop it from happening to begin with so these values don't skew the results, not just replace the negatives in the final output.

推荐答案

我的标准方法是将状态变量转换为不受限制的比例.对正变量执行此操作的最明显/标准方法是为log(x)而不是x的动力学写下方程.

My standard approach is to transform the state variables to an unconstrained scale. The most obvious/standard way to do this for positive variables is to write down equations for the dynamics of log(x) rather than of x.

例如，对于传染病流行的敏感传染病恢复(SIR)模型，当方程为dS/dt = -beta*S*I; dI/dt = beta*S*I-gamma*I; dR/dt = gamma*I时，我们会天真地将梯度函数写为

For example, with the Susceptible-Infected-Recovered (SIR) model for infectious disease epidemics, where the equations are dS/dt = -beta*S*I; dI/dt = beta*S*I-gamma*I; dR/dt = gamma*I we would naively write the gradient function as

gfun <- function(time, y, params) {
   g <- with(as.list(c(y,params)),
       c(-beta*S*I,
          beta*S*I-gamma*I,
          gamma*I)
       )
   return(list(g))
}

如果我们使log(I)而不是I作为状态变量(原则上我们也可以使用S来做到这一点，但是实际上S不太可能接近边界)，则我们有d(log(I))/dt = (dI/dt)/I = beta-gamma；其余等式需要使用exp(logI)来引用I.所以:

If we make log(I) rather than I be the state variable (in principle we could do this with S as well, but in practice S is much less likely to approach the boundary), then we have d(log(I))/dt = (dI/dt)/I = beta-gamma; the rest of the equations need to use exp(logI) to refer to I. So:

gfun_log <- function(time, y, params) {
   g <- with(as.list(c(y,params)),
       c(-beta*S*exp(logI),
          beta-gamma,
          gamma*exp(logI))
       )
   return(list(g))
}

(一次计算exp(logI)并存储/重用它比两次计算的效率更高……)

(it would be slightly more efficient to compute exp(logI) once and store/re-use it rather than computing it twice ...)

这篇关于用零替换模型(ODE系统)中的负值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！