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

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

问题是：在模拟期间，会有几次疫苗接种活动将人员转移到州内S陈述I。每个疫苗接种活动的特点是开始日期, an end date, its 覆盖率 and duration。我想做的是一旦time t落入区间开始日期 and end date对于一项疫苗接种活动，代码计算有效疫苗接种率（也m*n矩阵，基于覆盖率 and duration）并乘以S (m*n矩阵），得到过渡到状态的人员矩阵I。现在，我正在使用if()决定是否time t是介于一个开始日期 and a 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]
    }
  }
}

Here, irrate_matrix is the m*n有效疫苗接种率矩阵, m = 2是数量年龄组, n = 2是数量亚群, tbegin = c(5, 20, 35) is the 开始日期3 项疫苗接种活动中，tend = c(8, 23, 38) is the end date3 项疫苗接种活动中，covir = c(0.35, 0.25, 0.25, 0.225, 0.18, 0.13) is the 覆盖率每个亚群的每次疫苗接种（例如，covir[1] = 0.35 is the 覆盖率第一次疫苗接种的亚群1, while covir[4] = 0.225 is the 覆盖率第一次疫苗接种的亚群2) and 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)

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

times = seq(0, 50, 1)

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

我尝试过使用“事件”功能deSolvetpetzoldt 在这个问题中提到stackoverflow：根据时间更改参数值。然而，我发现改变一个参数矩阵并在每次有疫苗接种活动时都改变它对我来说很不方便。

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

• 如何改变我的irrate_matrix当有疫苗接种活动时将其设为非零矩阵，而在没有疫苗接种活动时将其设为零矩阵？ （必须为每次疫苗接种计算）

• 同时，如何通过避免随时迭代来让代码运行得更快tbegin[i] or 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)

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

实际上，CPU 时间已减少至 50%。相反，使用外部强迫approxTime而不是内在的if and 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解算器，但这对于您的完整模型可能有所不同。