## Sunday, December 13, 2015

### Gillespie, predator-prey simulation

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")``