r - 如何在 Gillespie 模拟中包含时变参数？

Question

我一直在模拟人口动态模型，并通过使其中一个参数的值随时间变化来添加一些环境随机性。

（更具体地说，我制作了将系统温度与系统内两个物种的生长速率相关联的热性能曲线。然后我随机抽取了一些温度以创建这些温度及其相应生长速率的向量。然后我设置模拟，使生长速率参数的值可以随系统温度随时间而变化。）

现在，我想使用 Gillespie 算法向我的系统添加一些人口统计随机性，特别是 R 中的 GillespieSSA 包。我在尝试将现有的环境随机实现与 ssa 函数采用的参数集成时遇到了麻烦。

这是我的环境随机实现：

  ri <- QuadEqn_1(temp = temp_sequence); rj <- QuadEqn_2(temp = temp_sequence) # Where QuadEqn_1 and _2 are the thermal performance curves that give the growth rate, when given the temperature of the system, "temp_sequence", which is a result of a random draw not shown here
  k <- 0.001; p <- 1 ; o <- 1000
  parms <- list(ri = ri, rj = rj, k = k, p = p, o = o)
  
  Antia_3sp_Model <- function(t,y,p1){
    tt <- floor(t) + 1
    Pi <- y[1]; Pj <- y[2]; I <- y[3]
    ri <- p1$ri[tt]; rj <- p1$rj[tt]; k <- p1$k; p <- p1$p; o <- p1$o # This is the line that allows the values of the growth rate parameters to change with time, tt
    
    dPi = ri*Pi - k*Pi*I # The model
    dPj = rj*Pj - k*Pj*I
    dI = p*I*(Pi/(Pi + o) + Pj/(Pj + o))
    list(c(dPi,dPj,dI))
  }
  N0 <- c(Pi = 1, Pj = 0, I = 1) # Initial values of the state variables
  TT <- round(seq(0.1, 50, 0.1), 1)
  eventdat <- data.frame(var = "Pj", time = 1, value = 1, method = "rep") # Allows one species to be introduced at different time points
  
  results <- lsoda(N0, TT, Antia_3sp_Model, p = parms, events=list(data=eventdat), verbose = TRUE)

使用 GillespieSSA 包需要一个倾向向量：

a <- c("ri*Pi",
           "k*Pi*I",
           "rj*Pj",
           "k*Pj*I",
           "p*I*(Pi/(Pi + o)",
           "p*I*(Pj/(Pj + o)")

还有一个状态变化矩阵：

nu <- matrix(c(+1, -1, 0, 0, 0, 0,
               0, 0, +1, -1, 0, 0,
               0, 0, 0, 0, +1, +1), nrow = 3, byrow = TRUE)

并且是使用 ssa 函数实现的，它看起来应该是这样的：

TestOutput <- ssa(x0, a, nu, parms1, TT, method, simName...)

所以，我想我的问题是：假设我通过倾向向量和状态变化矩阵将模型传递给 GillespieSSA，我如何才能像在环境随机实现中那样包含时变参数？

我对这个很迷茫，所以任何建议都将不胜感激:)

标签： rsimulationmodelingstochastic