ssa.exact {adaptivetau}

Exact stochastic simulation algorithm

Description

Simulates the trajectory of a continuous-time Markov process (aka. the Gillespie simulation algorithm).

Usage

ssa.exact(init.values, transitions, rateFunc, params, tf,
          reportTransitions = FALSE)

Arguments

init.values

Vector of initial values for all variables. This should be a simple one-dimensional vector of real numbers. IMPORTANT: variables must never take negative values.

transitions

Two-dimensional matrix of integers specifying how each state variable (rows) should be changed for a given transition (columns). Generally this will be a sparse matrix of primarily 1s and -1s.

rateFunc

R-function that returns instantaneous transition rates for each transition in the form a real-valued one-dimensional 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 transitions parameter above. This function must accept the following arguments:

• vector of current values for all state variables (in order used in the init.values argument above)

• parameters as supplied in argument to ssa.adaptivetau

• single real number giving the current time (all simulations start at t=0)

params

any R variable to be passed as-is to rateFunc, presumably specifying useful parameters.

tf

final time of simulation. Due to the approximate nature of the tau-leaping algorithm, the simulation may overshoot this time somewhat.

reportTransitions

default FALSE; whether to include a matrix of executed transitions per time interval in the output (alongside with the states at each time point).

Details

This function is supplied as a bonus with the adaptivetau package, since the C++ function that underlies this (exact) stochastic simulation algorithm is used in the (approximate) adaptive tau stochastic simulation as well.

The initial values, transition matrix, and transition rates must all be designed such that variables are always non-negative. The algorithm relies on this invariant.

Consider calling enableJIT(1) before running ssa.exact. In most circumstances, this should yield some speedup by byte-code compiling the rate function.

Value

Two-dimensional 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 the ‘reportTransitions’ option is used, then a list is returned with two elements. The first element is the ‘dynamics’ matrix described above. The second is a two-dimensional matrix called ‘transitions’ with a row for each timepoint and columns giving the number of times each transition was executed between the current time and the previous time.

Author(s)

Philip Johnson

See Also

This function is a bonus the comes along with the approximate (but faster) ssa.adaptivetau. For systems with sparse transition matrices, see helper function ssa.maketrans.

Examples

## Lotka-Volterra example
lvrates <- function(x, params, t) {
  with(params, {
    return(c(preygrowth*x["prey"],      ## prey growth rate
             x["prey"]*x["pred"]*eat,   ## prey death / predator growth rate
             x["pred"]*preddeath))      ## predator death rate
  })
}
params=list(preygrowth=10, eat=0.01, preddeath=10);
r=ssa.exact(c(prey = 1000, pred = 500),
            matrix(c(1,0, -2,1, 0,-1), nrow=2), lvrates, params, tf=2)
matplot(r[,"time"], r[,c("prey","pred")], type='l', xlab='Time', ylab='Counts')
legend("topleft", legend=c("prey", "predator"), lty=1:2, col=1:2)