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

Question

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

问题是：在模拟期间，会有几次疫苗接种活动，将状态S的人转移到状态I。每次疫苗接种 activity 都有一个 开始日期 、一个 结束日期 、 覆盖率 和 持续时间 。我想要做的是一旦 time t 落入 begin date 和 end date 的区间一次疫苗接种activity，代码计算有效疫苗接种率（也是一个m*n矩阵，基于覆盖率和持续时间)并用S乘以它(m*n matrix), 得到一个人的矩阵转移到状态I。现在，我正在使用 if() 来决定 time t 是否在 begin date 和 end date 之间:

#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 首次接种疫苗的 覆盖率，而 covir[4] = 0.225 是 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都接近[=28=的时间点] 或 tend[i]，模拟变得更慢，因为它在这个时间点迭代比任何其他时间点多得多的轮数。例如，一旦时间t来到tbegin[1] = 5，模型在时间点5迭代多轮。我附上了打印出这些迭代的屏幕截图 (screenshot1 and screenshot2)。我发现这就是为什么我的大模型现在需要很长时间运行ning 的原因。

我试过使用tpetzoldt在这个问题中提到的deSolve的“事件”功能。但是，我发现每次打疫苗都要改一个参数矩阵，对我来说很不方便activity.

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

如何在打疫苗的时候把我的irrate_matrix变成非零矩阵activity，在没有打疫苗的时候让它变成零矩阵？（每次接种都要计算）
同时，如何通过避免在任何tbegin[i]或tend[i]迭代多轮来使代码运行更快？（我想我不应该使用 if() 但我不知道我应该如何处理我的案例）
如果我需要使用“forcing”或“events”功能，你能告诉我如何在模型中有多个“forcing”/“events”吗？现在，我在更大的模型中使用了一个“事件”来将病毒引入人群，如：

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

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

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

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 求解器时，我观察到另一个相当大的加速，但这可能与您的完整模型不同。

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

R - deSolve package (ode function): change a matrix of parameters in SIR model according to time

r

ode

简化代码

基准