使用 ode 函数使 SEIR 模型参数取决于每个时间步的隔间发生率

Question

下面是我尝试使参数

r

取决于 SEIR 模型中

I

的发生率的代码。如果

(c(0, diff(I)))/10000

被替换为任何小于 1 的值，则模型将运行。但是，我希望参数

r

取决于

I

的值减去每个时间步上一个时间步的值

I

。

library(deSolve)

seirs <- function(t, x, parms)  {
  with(as.list(c(parms, x)), {
    dS= -lambda*S + rho*R
    dE= lambda*S - gamma_h*E
    dI=gamma_h*E - r*I
    dR=r*I - rho*R
    dCInc = lambda*S
    output <- c(dS, dE, dI, dR, dCInc)
    list(output)
  })
}

start<-c(S=5000, E=3000, I=2000,R=0, CInc = 0 )

parms <- c(lambda = 0.3, gamma_h = 1/20, r= (c(0, diff(I)))/10000, rho = 1/(1*365))
times <- seq(0, 10000, 1)
result <-ode(times=times, y=start, func=seirs,parms=parms)

Answer 1

参数通常是固定值。与时间相关的值可以通过多种方法实现：

1. 将依赖项集成到模型函数中，例如通过为参数添加一个额外的微分方程，
2. 使用强制函数并从外部“信号”进行插值，
3. 使用事件，修改状态变量。此方法可以与方法 #1 结合使用。

方法 #2 和 #3 在 deSolve 包的

?forcings

和

?events

帮助页面中进行了描述。

要使变量（或参数）依赖于当前时间步和前一个时间步的状态变量的差，首先需要定义微分方程系统中“时间步”的含义。根据定义，微分方程是连续的，因此理论上“时间步长”无限小。在数值模拟中，它有一个固定的小值，这取决于数值求解器。在大多数情况下，它会自动调整并且可能会有所不同。

因此，如果“时间步长”意味着“非常小”，那么 I 和前一个 I 之间的差值除以时间步长就是导数

dI

，请参见上面的方法#1。我认为这是首选方法

“时间步长”也可能意味着离散时间步长，例如1 天，假设

times

以天为单位。在这种情况下，可以使用延迟微分方程（DDE）。在 deSolve 中，可以使用

lagvalue

函数和

dede

求解器来实现。

以下代码包含此类 DDE 的示例。为了获得可信的图像，更改了参数值、状态变量和

times

向量。我还将差值从

I - I_tau

更改为

I_tau - I

以获得正值，并将

r

作为“辅助变量”添加到输出中。

请将此作为技术答案，因为我仍然不清楚原始问题的动机，即应该用这种方法对哪种生物过程进行建模。也许，上面的方法 #1 可能是更好的选择。

library(deSolve)

seirs <- function(t, x, parms)  {

  with(as.list(c(parms, x)), {
    if (t < tau)
      I_tau <- I # default value can also be set manually
    else
      I_tau <- lagvalue(t - tau, 1)

    r      <- k_r * (I_tau - I)
      
    dS     <- -lambda*S + rho*R
    dE     <- lambda*S - gamma_h*E
    dI     <- gamma_h*E - r*I
    dR     <- r*I - rho*R
    dCInc  <- lambda*S
    
    list(c(dS, dE, dI, dR, dCInc), r = r)
  })
}

start <- c(S=5000, E=3000, I=2000, R=0, CInc = 0)

parms <- c(lambda = 0.3, gamma_h = 1/20, k_r= 0.01, rho = 0.01, tau = 1)

times <- seq(0, 100, 1)
result <-dede(times=times, y=start, func=seirs, parms=parms)

plot(result)