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

如何解决如何在 Gillespie 模拟中包含时变参数？

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

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

现在，我想使用 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,+1,+1),nrow = 3,byrow = TRUE)

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

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

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

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