Approximate stochastic simulation algorithm
Description
Implements adaptive tauleaping approximation for simulating the trajectory of a continuoustime Markov process (see reference below).
Usage
1 2 3 4 5 
Arguments
init.values 
Vector of initial values for all variables. This should be a simple onedimensional vector of real numbers. IMPORTANT: variables must never take negative values. 
transitions 
One of two possible data types:
See the example below for details as well as

rateFunc 
R function that returns instantaneous transition rates for each
transition in the form a realvalued onedimensional vector with
length equal to the number of transitions. The order of these
rates must match the order in which transitions were specified in
the

params 
any R variable to be passed asis to rateFunc, presumably specifying useful parameters. 
tf 
final time of simulation. Due to the approximate nature of the tauleaping algorithm, the simulation may overshoot this time somewhat. 
jacobianFunc 
R function that returns Jacobian of transition rates. If supplied, enables the use of the implicit tauleaping algorithm (if appropriate; used in stiff systems). In contrast to typical definition of Jacobian (but consistent with the transition matrix parameter above), columns represent transitions (functions) and rows represent state variables (different partial derivative variables). jacobianFunc takes the same parameters as rateFunc. NOTE: the implicit tauleaping algorithm will reduce the observed variance for fastchanging state variables, although the mean will be correct. 
maxTauFunc 
R function that returns single real number giving the maximum tau leap allowable from the current state. Only executed if the adaptive tau algorithm wants to leap in a step greater than tl.params$maxtau (should only need to use if your rate functions have a dramatically discontinuous first derivative). Takes same parameters as rateFunc. 
deterministic 
Specify transitions to be treated as deterministic. If not NULL,
either a TRUE/FALSE vector of length equal to the number of
transitions or a integer vector of transition numbers to be
treated as deterministic (i.e. applying 
halting 
Specify transitions which, when executed, will cause the simulation
to halt. If not NULL, either a TRUE/FALSE vector of length equal to the
number of transitions or a integer vector of transition numbers to
be treated as halting (i.e. applying 
relratechange 
onedimensional vector of length equal to the number of state variables providing an upperbound to the ratio of amount that any transition rate will change given a corresponding change in the state variable. In other words, if variable i doubles, can we be assured that no transition will more than double (ratio of 1)? If not, then you need to set this variable to be greater than 1. 
tl.params 
List of various parameters to the tauleaping algorithm itself (best explained by reading original reference):

Details
The initial values, transition matrix & transition rates must all be designed such that variables are always nonnegative. The algorithm relies on this invariant.
See reference for details but, in brief the adaptive tauleaping algorithm dynamically switches between three methods for simulating events:
 exact
no approximation – executes a single transition at a time
 explicit tauleaping
subdivides transitions into those that might hit cause a variable to hit 0 (“critical”) and those that do not. Duration of time step picked dynamically with the goal of maximizing the step while minimized the change in the transition rates (the approximation assumes that these rates do not change). Noncritical transitions are advanced by a poissondistributed random variable; critical transitions are handled more like the exact algorithm.
 implicit tauleaping
in addition to dividing between critical and noncritical, also identifies transitions in quasiequilbrium (reversible pairs of transition that have roughly equal forwardbackward flow). Duration of time step picked on basis of noncritical, nonequilibrium transitions. This has the potential to greatly increase the timestep size for stiff systems. Similar idea to the explicit method but necesitates solving an implicit equation via Newton's method. Thus you must supply a function to calculate the Jacobian to enable this method.
All error messages that reference variables or transitions by number use 1based numbering where the first variable (or transition) is 1.
Consider calling enableJIT(1)
before running
ssa.adaptivetau. In most circumstances, this should yield some
speedup by bytecode compiling the rate function.
Value
If no halting transitions are specified, then a twodimensional matrix with rows for every timepoint at which the rateFunc was evaluated and columns giving the value for every state variable at that time. The first column specifies the time.
If halting transitions are specified, then a list with two elements. The first (‘dynamics’) is the same twodimensional matrix as without halting transitions. The second (‘haltingTransition’) gives the number of the transition that halted the simulation, or NA otherwise.
Note
Development of this package was supported in part by National Science Foundation award DBI0906041 and National Institute of Health award K99GM104158. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the NSF or NIH.
Author(s)
Philip Johnson
References
Cao Y, Gillespie DT, Petzold LR, The Journal of Chemical Physics, 2007
See Also
For systems with sparse transition matrices, see helper function
ssa.maketrans
. For a pureR implementation without C++
underneath, see GillespieSSA
. Also,
ssa.exact
exposes a R interface to the C++
implementation of the exact, nonapproximate simulation algorithm
(sometimes known as "Gillespie").
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41  ## Simple LotkaVolterra example
# We have three potential transitions:
transitions = list(c(prey = +1), # prey grow
c(prey = 2, pred = +1), # predation
c(pred = 1)) # predator dies
# Function to calculate transition rates, given variables and parameters
lvrates < function(x, params, t) {
return(c(params$preygrowth*x["prey"], # rate of prey growing
x["prey"]*x["pred"]*params$eat, # rate of predation
x["pred"]*params$preddeath)) # rate of predators dying
}
# Set the LotkaVolterra parameters
params = list(preygrowth=10, eat=0.01, preddeath=10);
# Set the random seed (only necessary if you want to reproduce results)
set.seed(4)
# Perform the stochastic simulation!
r=ssa.adaptivetau(c(prey = 1000, pred = 500),
transitions, lvrates, params, tf=12)
# Plot the results
matplot(r[,"time"], r[,c("prey","pred")], type='l', xlab='Time',
ylab='Counts (log scale)', log='y')
legend("bottomleft", legend=c("prey", "predator"), lty=1:2, col=1:2)
# However, if you are interested in very finescale variance, perhaps the
# default parameters smooth too much. Try reducing the tl.param epsilon
# to see a better approximation:
s=ssa.adaptivetau(c(prey = 1000, pred = 500),
transitions, lvrates, params, tf=12,
tl.params = list(epsilon=0.01)) # reduce "epsilon"
par(mfrow=c(2,1));
matplot(r[r[,"time"]<2,"time"], r[r[,"time"]<2,c("prey","pred")],
type='l', xlab='Time', ylab='Counts', main="Original (epsilon=default)")
matplot(s[s[,"time"]<2,"time"], s[s[,"time"]<2,c("prey","pred")],
type='l', xlab='Time', ylab='Counts', main="Finescale (epsilon=0.01)")
