Блог

Я моделировал модель динамики численности населения и добавил некоторую стохастичность окружающей среды, изменив значение одного из параметров во времени.

(Чтобы быть более конкретным, я сделал кривые тепловых характеристик, которые связывают температуру системы со скоростью роста двух видов в системе. Затем я произвольно произвел выборку некоторых температур, чтобы создать вектор этих температур и соответствующих им темпов роста. Затем я настроил моделирование таким образом, чтобы значение параметров скорости роста могло изменяться с течением времени в зависимости от температуры системы.)

Теперь я хочу добавить некоторую демографическую стохастичность в свою систему, используя алгоритм Гиллеспи, и, в частности, пакет GillespieSSA в R. Я столкнулся с проблемами, пытаясь интегрировать мою существующую реализацию стохастики окружающей среды с аргументами, которые принимают функции 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...)
 

Итак, я предполагаю, что мой вопрос таков: учитывая, что я передаю модель Гиллеспи через вектор склонности и матрицу изменения состояния, как я могу включить изменяющийся во времени параметр, как это было в моей экологически стохастической реализации?

Я довольно запутался в этом вопросе, так что любые предложения будут весьма признательны 🙂