R - deSolve 包（ode 函数）：根据时间更改 SIR 模型中的参数矩阵

Question

我正在尝试使用 deSolve 包中的函数 ode 来模拟病毒在人群中的传播。我的模型的基础是 SIR 模型，我在这里发布了一个更简单的模型演示，它仅包含三个状态 S（易感）、I（传染性）和 R（已恢复）。在我的代码中，每个州都由 m*n 矩阵 表示，因为我的人口中有 m 个年龄组 和 n 个子群体。 p>

问题是：在模拟期间，会有几次疫苗接种活动将状态S的人转移到状态I。每个疫苗接种活动的特点是开始日期、结束日期、覆盖率和持续时间。我想要做的是一旦 时间 t 落入一个疫苗接种活动的 开始日期 和 结束日期 的间隔，代码计算有效疫苗接种率（也是一个m*n矩阵，基于覆盖率和持续时间）并将其与S（m*n矩阵），得到一个人的矩阵转换到状态I。现在，我正在使用if() 来确定时间 t 是否介于 开始日期 和 结束日期 之间：

#initialize the matrix of effective vaccination rate 
irrate_matrix = matrix(data = rep(0, m*n), nrow = m, ncol = n) 

for (i in 1:length(tbegin)){
  if (t>=tbegin[i] & t<=tend[i]){
    for (j in 1:n){
      irrate_matrix[, j] = -log(1-covir[(j-1)*length(tbegin)+i])/duration[i]
    }
  }
}

这里，irrate_matrix是m*n有效接种率矩阵，m = 2是年龄组数，n = 2是亚群，tbegin = c(5, 20, 35) 是 3 次疫苗接种活动的开始日期，tend = c(8, 23, 38) 是 3 次疫苗接种活动的结束日期，covir = c(0.35, 0.25, 0.25, 0.225, 0.18, 0.13) 是每个亚群的每次疫苗接种的覆盖率（例如，covir[1] = 0.35 是 subpopulation1 的第一次疫苗接种的覆盖率，而@987654335 @ 是 subpopulation2 的第一次疫苗接种的覆盖率，duration = c(4, 4, 4) 是每次疫苗接种的持续时间（以天为单位）。

在计算irrate_matrix 后，我将其转化为导数，因此我有：

dS = as.matrix(b*N) - as.matrix(irrate_matrix*S) - as.matrix(mu*S)
dI = as.matrix(irrate_matrix*S) - as.matrix(gammaS*I) - as.matrix(mu*I)
dR = as.matrix(gammaS*I) - as.matrix(mu*R)

我想从第 0 天到第 50 天进行模拟，以 1 天为单位，因此：

times = seq(0, 50, 1)

我的代码的当前问题是：每次 时间 t 到达接近 tbegin[i] 或 tend[i] 的时间点时，模拟变得更多较慢，因为它在这个时间点迭代比在任何其他时间点更多的轮次。例如，一旦 时间 t 达到tbegin[1] = 5，模型会在时间点 5 迭代许多轮。我附上了打印这些迭代的屏幕截图（screenshot1 和screenshot2）。我发现这就是为什么我的更大模型现在需要很长时间运行的原因。

我已经尝试使用 tpetzoldt 在这个问题stackoverflow: change the value of a parameter as a function of time 中提到的deSolve 的“事件”功能。但是，我发现每次有疫苗接种活动时都更改参数矩阵并更改它对我来说很不方便。

我正在寻找以下方面的解决方案：

• 如何在有疫苗接种活动时将我的irrate_matrix 更改为非零矩阵，并在没有疫苗接种时让它为零矩阵？ （每次接种都要计算）

• 同时，如何避免在任何tbegin[i] 或tend[i] 处进行多次迭代，从而使代码运行得更快？ （我认为我不应该使用if()，但我不知道我应该如何处理我的案例）

• 如果我需要使用“强制”或“事件”功能，您能否告诉我如何在模型中设置多个“强制”/“事件”？现在，我在更大的模型中使用了一个“事件”来向人群引入病毒，例如：

virusevents = data.frame(var = "I1", time = 2, value = 1, method = "add")

欢迎任何好主意，非常感谢直接提供一些代码！提前谢谢！

作为参考，我在这里发布整个演示：

library(deSolve)

##################################
###(1) define the sir function####
##################################

sir_basic <- function (t, x, params) 
{ # retrieve initial states
  S = matrix(data = x[(0*m*n+1):(1*m*n)], nrow = m, ncol = n)
  I = matrix(data = x[(1*m*n+1):(2*m*n)], nrow = m, ncol = n)
  R = matrix(data = x[(2*m*n+1):(3*m*n)], nrow = m, ncol = n)
  
  with(as.list(params), {
    N = as.matrix(S + I + R) 
    
    # print out current iteration
    print(paste0("Total population at time ", t, " is ", sum(N)))
    
    # calculate irrate_matrix by checking time t
    irrate_matrix = matrix(data = rep(0, m*n), nrow = m, ncol = n)
    for (i in 1:length(tbegin)){
      if (t>=tbegin[i] & t<=tend[i]){
        for (j in 1:n){
          irrate_matrix[, j] = -log(1-covir[(j-1)*length(tbegin)+i])/duration[i]
        }
      }
    }
    
    # derivatives
    dS = as.matrix(b*N) - as.matrix(irrate_matrix*S) - as.matrix(mu*S)
    dI = as.matrix(irrate_matrix*S) - as.matrix(gammaS*I) - as.matrix(mu*I)
    dR = as.matrix(gammaS*I) - as.matrix(mu*R)
    derivatives <- c(dS, dI, dR)
    list(derivatives)
  })
}

##################################
###(2) characterize parameters####
##################################

m = 2 # the number of age groups
n = 2 # the number of sub-populations

tbegin = c(5, 20, 35) # begin dates
tend = c(8, 23, 38) # end dates
duration = c(4, 4, 4) # duration
covir = c(0.35, 0.25, 0.25, 0.225, 0.18, 0.13) # coverage rates

b = 0.0006 # daily birth rate
mu = 0.0006 # daily death rate
gammaS = 0.05 # transition rate from I to R

parameters = c(m = m, n = n,
               tbegin = tbegin, tend = tend, duration = duration, covir = covir,
               b = b, mu = mu, gammaS = gammaS)

##################################
#######(3) initial states ########
##################################

inits = c(
  S = c(20000, 40000, 10000, 20000),
  I = rep(0, m*n),
  R = rep(0, m*n)
)

##################################
#######(4) run simulations########
##################################

times = seq(0, 50, 1)

traj <- ode(func  = sir_basic, 
            y = inits,
            parms = parameters,
            times = times)

plot(traj)

Answer 1

矩阵和向量的元素明智操作相同，因此as.matrix 转换是多余的，因为没有使用true 矩阵乘法。与rep 相同：无论如何都会回收零。

实际上，CPU 时间已经减少到 50%。相比之下，使用带有approxTime 的外部强制而不是内部if 和for 会使模型变慢（未显示）。

简化代码

sir_basic2 <- function (t, x, params)
{ # retrieve initial states
  S = x[(0*m*n+1):(1*m*n)]
  I = x[(1*m*n+1):(2*m*n)]
  R = x[(2*m*n+1):(3*m*n)]

  with(as.list(params), {
    N = S + I + R

    # print out current iteration
    #print(paste0("Total population at time ", t, " is ", sum(N)))

    # calculate irrate_matrix by checking time t
    irrate_matrix = matrix(data = 0, nrow = m, ncol = n)
    for (i in 1:length(tbegin)){
      if (t >= tbegin[i] & t <= tend[i]){
        for (j in 1:n){
          irrate_matrix[, j] = -log(1-covir[(j-1) * length(tbegin)+i])/duration[i]
        }
      }
    }

    # derivatives
    dS = b*N - irrate_matrix*S - mu*S
    dI = irrate_matrix*S - gammaS*I - mu*I
    dR = gammaS*I - mu*R
    list(c(dS, dI, dR))
  })
}

基准测试

每个模型版本运行 10 次。模型sir_basic 是原始实现，其中print 行被禁用以进行公平比较。

system.time(
  for(i in 1:10)
    traj <- ode(func  = sir_basic,
            y = inits,
            parms = parameters,
            times = times)
)

system.time(
  for(i in 1:10)
    traj2 <- ode(func  = sir_basic2,
            y = inits,
            parms = parameters,
            times = times)
)

plot(traj, traj2)
summary(traj - traj2)

当我使用 method="adams" 而不是默认的 lsoda 求解器时，我观察到另一个显着的加速，但这对于您的完整模型可能会有所不同。