This function creates a function for advancing the state of an SPN model using the Gillespie algorithm. The resulting function (closure) can be used in conjunction with other functions (such as simTs) for simulating realisations of SPN models.
StepGillespie(N)An R function which can be used to advance the state of the SPN model N by using the Gillespie algorithm. The function closure has interface function(x0,t0,deltat,...), where x0 and t0 represent the initial state and time, and deltat represents the amount of time by which the process should be advanced. The function closure returns a vector representing the simulated state of the system at the new time.
An R list with named components representing a stochastic
Petri net (SPN). Should contain N$Pre, a matrix representing
the LHS stoichiometries, N$Post, a matrix representing the
RHS stoichiometries, and N$h, a function representing the
rates of the reaction processes. N$h should have
first argument x, a vector representing the
current state of the system, and
second
argument t, a scalar representing the current simulation time
(in the typical time-homogeneous case, N$h will ignore this
argument).
N$h may possess additional arguments, representing reaction rates, for example. N does not need to contain an initial marking, N$M. N$M will be ignored by most functions which use the resulting function closure.
StepEulerSPN, StepGillespie1D,
simTs, simTimes, simSample, StepFRM,
StepPTS, StepCLE
# load up the Lotka-Volterra (LV) model
data(spnModels)
LV
# create a stepping function
stepLV = StepGillespie(LV)
# step the function
print(stepLV(c(x1=50,x2=100),0,1))
# simulate a realisation of the process and plot it
out = simTs(c(x1=50,x2=100),0,100,0.1,stepLV)
plot(out)
plot(out,plot.type="single",lty=1:2)
# simulate a realisation using simTimes
times = seq(0,100,by=0.1)
plot(ts(simTimes(c(x1=50,x2=100),0,times,stepLV),start=0,deltat=0.1),plot.type="single",lty=1:2)
# simulate a realisation at irregular times
times = c(0,10,20,50,100)
out2 = simTimes(c(x1=50,x2=100),0,times,stepLV)
print(out2)
Run the code above in your browser using DataLab