How to include time-varying parameters in a Gillespie simulation?

Question

I've been simulating a population dynamics model, and added in some environmental stochasticity by making the value of one of the parameters time-varying.

(To be more specific, I made thermal performance curves that relate the temperature of the system to the growth rates of two species within the system. I then randomly sampled some temperatures to create a vector of those temperatures and their corresponding growth rates. I then set up the simulation such that the value of the growth rate parameters could change with the temperature of the system over time.)

Now, I want to add some demographic stochasticity to my system using the Gillespie algorithm, and specifically, the GillespieSSA package in R. I've run into trouble trying to integrate my existing environmentally stochastic implementation with the arguments the ssa functions takes.

This is my environmentally stochastic implementation:

  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)

Using the GillespieSSA package requires a propensity vector:

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

And a state change matrix:

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)

And is implemented using the ssa function, which should look something like this:

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

So, I guess my question is: given that I pass GillespieSSA the model through the propensity vector and the state change matrix, how can I include a time-varying parameter, as I have in my environmentally stochastic implementation?

I'm pretty lost on this one, so any suggestions would be greatly appreciated :)