Gillespie, predator-prey simulation

GillespieSSA test. Predator-Prey

H Qin

December 13, 2015

Example from Pineda-Krch08 J of Statistical Software.

rm(list=ls())
library(GillespieSSA)
parms = c(b=2, d=1, K=1000, alpha=0.007, w=0.0035, c=2,g=2)
x0 = c(N=1000, P=100)
nu = matrix(c(+1,-1,-1, 0,0,
              0,  0, 0,+1,-1), nrow=2, byrow=TRUE)
a = c("b*N", 
      "(d+(b-d)*N/K)*N", 
      "alpha/(1+w*N)*N*P",
      "c*alpha/(1+w*N)*N*P",
      "g*P")
tf = 100
method = "D"
simName = "Predator-prey model"

parms

##       b       d       K   alpha       w       c       g 
## 2.0e+00 1.0e+00 1.0e+03 7.0e-03 3.5e-03 2.0e+00 2.0e+00

nu 

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1   -1   -1    0    0
## [2,]    0    0    0    1   -1

a

## [1] "b*N"                 "(d+(b-d)*N/K)*N"     "alpha/(1+w*N)*N*P"  
## [4] "c*alpha/(1+w*N)*N*P" "g*P"

out = ssa(x0, a, nu, parms, tf, method, simName, 
          verbose = TRUE, consoleInterval = 10, maxWallTime = 30)

## Running D method with console output every 10 time step
## Start wall time: 2015-12-13 20:44:52...
## t=0 : 1000,100
## (4.082s) t=10 : 639,313
## (10.931s) t=20.00066 : 123,499
## (18.091s) t=30.00005 : 112,193
## (26.468s) t=40.00018 : 271,95
## t=43.15012 : 374,361
## tf: 43.15012
## TerminationStatus: maxWallTime
## Duration: 30.019 seconds
## Method: D
## Nr of steps: 96521
## Mean step size: 0.0004470542+/-0.0006694862
## End wall time: 2015-12-13 20:45:22
## --------------------

summary(out$data)

##        V1               N                P      
##  Min.   : 0.000   Min.   :  46.0   Min.   :  9  
##  1st Qu.: 9.575   1st Qu.: 294.0   1st Qu.: 92  
##  Median :19.063   Median : 496.0   Median :206  
##  Mean   :20.357   Mean   : 487.3   Mean   :237  
##  3rd Qu.:31.002   3rd Qu.: 677.0   3rd Qu.:361  
##  Max.   :43.150   Max.   :1014.0   Max.   :622

plot( out$data[,2] ~ out$data[,1], ylab = "N")

plot( out$data[,3] ~ out$data[,1], ylab = "Predator")