Delayed stochastic algorithm extend the Gillespie stochastic simulation algorithm (SSA) to a system with delays. A delayed reaction consists of an exponential waiting time as initiation step with a subsequent delay time. We implemented both the Delayed Direct Method (Cai 2007, J. Chem. Phys) and the Delayed Next Reaction Method (Anderson 2007, J. Chem. Phys)
"""
Written by TR Maarleveld, Amsterdam, The Netherlands
E-mail: tmd200@users.sourceforge.net
Last Change: October 08, 2014
"""

import stochpy
smod = stochpy.SSA()
# required for iPython Notebook inline plotting
%matplotlib inline 

We start with a SSA without delay and subsequently do two types of delayed SSAss: a "consuming" delay" followed by a "nonconsuming" delay. For consuming delayed reactions, the reactants are consumed at initiation and for nonconsuming delayed reactions, the reactants are consumed at completion. In both cases, the products are updated at completion. By default, delayed reactions are “consuming delayed reactions”. 

In this example, we set a fixed delay of 5 time units for reaction R1. Note that we use the delayed direct method, but that we can also use the delayed next reaction method.
smod.Model('Isomerization.psc')
smod.DoStochSim(mode='time',end=10,trajectories=1000)
smod.GetRegularGrid()
smod.PlotAverageSpeciesTimeSeries()

smod.SetDelayParameters({'R1':('fixed',5)})
smod.DoDelayedStochSim(mode='time',end=10,trajectories=1000)
smod.GetRegularGrid()
smod.PlotAverageSpeciesTimeSeries()

smod.SetDelayParameters({'R1':('fixed',5)},nonconsuming_reactions=['R1'])
smod.DoDelayedStochSim(mode='time',end=10,trajectories=1000)
smod.GetRegularGrid()
smod.PlotAverageSpeciesTimeSeries()
We now use the "CellDivision" to explore delayed SSAs in a bit more detail. For example, we here use three different types of delays (fixed, lognormal, and gamma) and mix both consuming (R3, R4) and nonconsuming delays (R6).

Stochastic simulations are done (i) with the delayed direct method and (ii) with the delayed next reaction method. The results will not be identical, but this is normal for stochastic simulations.
smod.Model('CellDivision.psc')
smod.SetDelayParameters({'R3':('fixed',2), 'R4':('lognormal',3,1), 'R6':('gamma',4,1)}, nonconsuming_reactions = ['R6'])
# Same result when using index instead of reaction name (with 'R1' = index 0).
smod.SetDelayParameters({2:('fixed',2),3:('lognormal',3,1),5:('gamma',4,1)}, nonconsuming_reactions = [5])
smod.DoDelayedStochSim()
smod.PlotSpeciesTimeSeries()
smod.DoDelayedStochSim(method = 'delayedNRM')   # Delay parameters are remembered
smod.PlotSpeciesTimeSeries()
Note that the correct high-level function must be used to perform delayed stochastic simulations. We illustrate this below.
# DoStochSim will switch automatically back to normal direct method if delayed method given.
smod.DoStochSim(method = 'delayeddirect') 
smod.PlotSpeciesTimeSeries()
# Similarly, DoDelayedStochSim will switch automatically back to direct delayed method if normal method given.
smod.DoDelayedStochSim(method = 'nrm')  
smod.PlotSpeciesTimeSeries()