ssa  R Documentation 
Main interface function to the implemented SSA methods. Runs a single realization of a predefined system.
ssa( x0, # initial state vector a, # propensity vector nu, # statechange matrix parms = NULL, # model parameters tf, # final time method = ssa.d(), # SSA method simName = "", tau = 0.3, # deprecated f = 10, # deprecated epsilon = 0.03, # deprecated nc = 10, # deprecated hor = NA_real_, # deprecated dtf = 10, # deprecated nd = 100, # deprecated ignoreNegativeState = TRUE, consoleInterval = 0, censusInterval = 0, verbose = FALSE, maxWallTime = Inf )
x0 
numerical vector of initial states where the component elements
must be named using the same notation as the corresponding state variable in
the propensity vector, 
a 
character vector of propensity functions where state variables
correspond to the names of the elements in 
nu 
numerical matrix of change if the number of individuals in each state (rows) caused by a single reaction of any given type (columns). 
parms 
named vector of model parameters. 
tf 
final time. 
method 
an SSA method, the valid options are:

simName 
optional text string providing an arbitrary name/label for the simulation. 
tau 
[DEPRECATED], see 
f 
[DEPRECATED], see 
epsilon 
[DEPRECATED], see 
nc 
[DEPRECATED], see 
hor 
[DEPRECATED], see 
dtf 
[DEPRECATED], see 
nd 
[DEPRECATED], see 
ignoreNegativeState 
boolean object indicating if negative state
values should be ignored (this can occur in the 
consoleInterval 
(approximate) interval at which 
censusInterval 
(approximate) interval between recording the state of the system.
If 
verbose 
boolean object indicating if the status of the simulation
simulation should be displayed on the console. If 
maxWallTime 
maximum wall time duration (in seconds) that the simulation is allowed to run for before terminated. This option is useful, in particular, for systems that can end up growing uncontrolably. 
Although ssa
can be invoked by only specifying the system
arguments (initial state vector x0
, propensity vector a
,
statechange matrix nu
), the final time (tf
), and the
SSA method to use, substantial improvements in speed and accuracy
can be obtained by adjusting the additional (and optional) ssa
arguments. By default ssa
(tries to) use conservative default values
for the these arguments, prioritizing computational accuracy over
computational speed. These default values are, however, not fool
proof for the approximate methods, and occasionally one will have to hand
tweak them in order for a stochastic model to run appropriately.
Returns a list object with the following elements,
data
: a numerical matrix object of the simulation time series where the first column is the time vector and subsequent columns are the state frequencies.
stats
: sublist object with elements containing various simulation statistics. The of the sublist are:
stats$startWallTime
: start wall clock time (YYYYmmdd HH:MM:SS).
stats$endWallTime
: end wall clock time (YYYYmmdd HH:MM:SS).
stats$elapsedWallTime
: elapsed wall time in seconds.
stats$terminationStatus
: string vector listing the reason(s) for the
termination of the realization in 'plain words'. The possible termination statuses are:
finalTime
= if the simulation reached the maximum simulation time tf
,
extinction
= if the population size of all states is zero,
negativeState
= if one or several states have a negative population size (can occur in the ETL method),
zeroProp
= if all the states have a zero propensity function,
maxWallTime
= if the maximum wall time has been reached. Note the termination status may have more than one message.
'stats$nSteps“ total number of time steps (or tauleaps) executed.
stats$meanStepSize
: mean step (or tauleap) size.
stats$sdStepSize
: one standard deviation of the step (or tauleap) size.
stats$SuspendedTauLeaps
: number of steps performed using the Direct method due to OTL
suspension (only applicable for the OTL
method).
arg$...
: sublist with elements containing all the arguments and their values used to invoke ssa
(see Usage and Arguments list above).
In order to invoke SSA the stochastic
model needs at least four components, the initial state vector (x0
),
statechange matrix (nu
), propensity vector (a
), and the final
time of the simulation (tf
). The initial state vector defines the
population sizes in all the states at t=0, e.g. for a system with two
species X1
and X2
where both have an initial population size
of 1000 the initial state vector is defined as x0 < c(X1=1000,X2=1000)
. The elements of the vector have to be labelled using
the same notation as the state variables used in the propensity functions.
The statechange matrix defines the change in the number of individuals in
each state (rows) as caused by one reaction of a given type (columns). For
example, the statechange matrix for system with the species S1
and S2 with two reactions
S1 c1> S2
S2 c2> 0
is defined as nu < matrix(c(1,0,+1,1),nrow=2,byrow=TRUE)
where
c1 and c2 are the per capita reaction probabilities.
The propensity vector, a
, defines the probabilities that a particular
reaction will occur over the next infinitesimal time interval [t,t+dt]. For example, in the previous example the
propensity vector is defined as a < c("c1*X1","c2*X2")
. The
propensity vector consists of character elements of each reaction's
propensity function where each state variable requires the corresponding
named element label in the initial state vector (x0
).
Irreversible isomerization: Perhaps the simplest model that can be formulated using the SSA is the irreversible isomerization (or radioactive decay) model. This model is often used as a first pedagogic example to illustrate the SSA (see e.g. Gillespie 1977). The deterministic formulation of this model is
dX/dt=cX
where the single reaction channel is
S c> 0.
By setting X_0=1000 and c=0.5 it is now simple to define this model and run it for 10 time steps using the Direct method,
out < ssa(x0=c(X=1000),a=c("c*X"),nu=matrix(1),parms=c(c=0.5),tf=10)
The resulting time series can then be displayed by,
ssa.plot(out)
Selecting the appropriate SSA method is a tradeoff between computational speed, accuracy of the results, and which SSA actually works for a given scenario. This depends on the characteristics of the defined system (e.g. number of reaction channels, number of species, and the absolute and relative magnitude of the propensity functions). Not all methods are appropriate for all models. When selecting a SSA method all of these factors have to be taken into consideration. The various tauleap methods accept a number of additional arguments. While the default values of these arguments may work for some scenarios they may have to be adjusted for others. The default values for the tauleap methods are conservative in terms of computational speed and substantial increase in efficiency may be gained by optimizing their values for a specific system.
GillespieSSApackage, ssa.d()
, ssa.etl()
, ssa.btl()
, ssa.otl()
## Irreversible isomerization ## Large initial population size (X=1000) parms < c(c=0.5) x0 < c(X=10000) a < c("c*X") nu < matrix(1) out < ssa(x0,a,nu,parms,tf=10,method=ssa.d(),simName="Irreversible isomerization") # Direct method plot(out$data[,1],out$data[,2]/10000,col="red",cex=0.5,pch=19) ## Smaller initial population size (X=100) x0 < c(X=100) out < ssa(x0,a,nu,parms,tf=10,method=ssa.d()) # Direct method points(out$data[,1],out$data[,2]/100,col="green",cex=0.5,pch=19) ## Small initial population size (X=10) x0 < c(X=10) out < ssa(x0,a,nu,parms,tf=10,method=ssa.d()) # Direct method points(out$data[,1],out$data[,2]/10,col="blue",cex=0.5,pch=19) ## Logistic growth parms < c(b=2, d=1, K=1000) x0 < c(N=500) a < c("b*N", "(d+(bd)*N/K)*N") nu < matrix(c(+1,1),ncol=2) out < ssa(x0,a,nu,parms,tf=10,method=ssa.d(),maxWallTime=5,simName="Logistic growth") ssa.plot(out) ## KermackMcKendrick SIR model parms < c(beta=0.001, gamma=0.1) x0 < c(S=499,I=1,R=0) a < c("beta*S*I","gamma*I") nu < matrix(c(1,0,+1,1,0,+1),nrow=3,byrow=TRUE) out < ssa(x0,a,nu,parms,tf=100,method=ssa.d(),simName="SIR model") ssa.plot(out) ## Lotka predatorprey model parms < c(c1=10, c2=.01, c3=10) x0 < c(Y1=1000,Y2=1000) a < c("c1*Y1","c2*Y1*Y2","c3*Y2") nu < matrix(c(+1,1,0,0,+1,1),nrow=2,byrow=TRUE) out < ssa(x0,a,nu,parms,tf=100,method=ssa.etl(),simName="Lotka predatorprey model") ssa.plot(out)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.