GillespieSSA (version 0.6.2)

ssa: Invoking the stochastic simulation algorithm

Description

Main interface function to the implemented SSA methods. Runs a single realization of a predefined system.

Usage

ssa(
                   x0,            # initial state vector
                    a,            # propensity vector
                   nu,            # state-change 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
)

Arguments

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.

a

character vector of propensity functions where state variables correspond to the names of the elements in x0.

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:

• ssa.d() --- Direct method (default method),

• ssa.etl() - Explicit tau-leap,

• ssa.btl() --- Binomial tau-leap, or

• ssa.otl() --- Optimized tau-leap.

simName

optional text string providing an arbitrary name/label for the simulation.

tau

[DEPRECATED], see ssa.etl()

f

[DEPRECATED], see ssa.btl()

epsilon

[DEPRECATED], see ssa.otl()

nc

[DEPRECATED], see ssa.otl()

hor

[DEPRECATED], see ssa.otl()

dtf

[DEPRECATED], see ssa.otl()

nd

[DEPRECATED], see ssa.otl()

ignoreNegativeState

boolean object indicating if negative state values should be ignored (this can occur in the etl method). If ignoreNegativeState=TRUE the simulation finishes gracefully when encountering a negative population size (i.e. does not throw an error). If ignoreNegativeState=FALSE the simulation stops with an error message when encountering a negative population size.

consoleInterval

(approximate) interval at which ssa produces simulation status output on the console (assumes verbose=TRUE). If consoleInterval=0 console output is generated each time step (or tau-leap). If consoleInterval=Inf no console output is generated. Note, verbose=FALSE disables all console output. Console output drastically slows down simulations.

censusInterval

(approximate) interval between recording the state of the system. If censusInterval=0 \((t,x)\) is recorded at each time step (or tau-leap). If censusInterval=Inf only \((t_0,x_0)\) and \((t_f,x_t)\) is recorded. Note, the size of the time step (or tau-leaps) ultimately limits the interval between subsequent recordings of the system state since the state of the system cannot be recorded at a finer time interval the size of the time steps (or tau-leaps).

verbose

boolean object indicating if the status of the simulation simulation should be displayed on the console. If verbose=TRUE the elapsed wall time and \((t,x)\) is displayed on the console every consoleInterval time step and a brief summary is displayed at the end of the simulation. If verbose=FALSE the simulation runs entirely silent (overriding consoleInterval). Verbose runs drastically slows down simulations.

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.

Value

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: sub-list object with elements containing various simulation statistics. The of the sub-list are:

• stats$startWallTime: start wall clock time (YYYY-mm-dd HH:MM:SS).

• stats$endWallTime: end wall clock time (YYYY-mm-dd 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 tau-leaps) executed.

• stats$meanStepSize: mean step (or tau-leap) size.

• stats$sdStepSize: one standard deviation of the step (or tau-leap) size.

• stats$SuspendedTauLeaps: number of steps performed using the Direct method due to OTL suspension (only applicable for the OTL method).

• arg$...: sub-list with elements containing all the arguments and their values used to invoke ssa (see Usage and Arguments list above).

Preparing a run

In order to invoke SSA the stochastic model needs at least four components, the initial state vector (x0), state-change 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 state-change 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 state-change matrix for system with the species \(S_1\) and \(S_2\) with two reactions $$S_1 \stackrel{c_1}{\longrightarrow} S_2$$ $$S_2 \stackrel{c_2}{\longrightarrow} 0$$

is defined as nu <- matrix(c(-1,0,+1,-1),nrow=2,byrow=TRUE) where \(c_1\) and \(c_2\) 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 \(\left[ t,t+dt \right]\). 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).

Example

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

$$\frac{dX}{dt}=-cX$$

where the single reaction channel is

$$S \stackrel{c}{\longrightarrow} 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)

Details

Although ssa can be invoked by only specifying the system arguments (initial state vector x0, propensity vector a, state-change 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.

See Also

GillespieSSA-package, ssa.d(), ssa.etl(), ssa.btl(), ssa.otl()

Examples

# NOT RUN {
## 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+(b-d)*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)

## Kermack-McKendrick 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 predator-prey 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 predator-prey model")
ssa.plot(out)

# }