StochPy's delayed stochastic simulation algorithms (delayed SSAs)

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)
In [1]:
"""
Written by TR Maarleveld, Amsterdam, The Netherlands
E-mail: [email protected]
Last Change: October 08, 2014
"""

import stochpy
smod = stochpy.SSA()
# required for iPython Notebook inline plotting
%matplotlib inline 
Info: Direct method is selected to perform stochastic simulations.
Parsing file: /home/timo/Stochpy/pscmodels/ImmigrationDeath.psc
Info: No reagents have been fixed

#######################################################################
#                                                                     #
#            Welcome to the interactive StochPy environment           #
#                                                                     #
#######################################################################
#  StochPy: Stochastic modeling in Python                             #
#  http://stochpy.sourceforge.net                                     #
#  Copyright(C) T.R Maarleveld, B.G. Olivier, F.J Bruggeman 2010-2014 #
#  DOI: 10.1371/journal.pone.0079345                                  #
#  Email: [email protected]                                #
#  VU University, Amsterdam, Netherlands                              #
#  Centrum Wiskunde Informatica, Amsterdam, Netherlands               #
#  StochPy is distributed under the BSD licence.                      #
#######################################################################

Version 2.2.0
Output Directory: /home/timo/Stochpy
Model Directory: /home/timo/Stochpy/pscmodels
Warning: Figure freezing? Try a different matplotlib backend (stochpy.plt.switch_backend) and/or set IsInteractive to False (see user guide)
Info: Direct method is selected to perform stochastic simulations.
Parsing file: /home/timo/Stochpy/pscmodels/ImmigrationDeath.psc
Info: No reagents have been fixed

Example 1:

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.
In [2]:
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()
Parsing file: /home/timo/Stochpy/pscmodels/Isomerization.psc
Info: No reagents have been fixed
Info: 1000 trajectories are generated
Info: Time simulation output of the trajectories is stored at Isomerization(trajectory).dat in directory: /home/timo/Stochpy/temp
Info: Simulation time: 6.27937412262                                                   
*** WARNING ***: an invalid method (Direct) was selected. Switching to the Delayed Direct Method.
Info: Delayed Direct Method is selected to perform delayed stochastic simulations.
Parsing file: /home/timo/Stochpy/pscmodels/Isomerization.psc
Info: No reagents have been fixed
Info: 1000 trajectories are generated
Info: Time simulation output of the trajectories is stored at Isomerization(trajectory).dat in directory: /home/timo/Stochpy/temp
Info: Simulation time: 9.16363191605                                                   
Info: 1000 trajectories are generated
Info: Time simulation output of the trajectories is stored at Isomerization(trajectory).dat in directory: /home/timo/Stochpy/temp
Info: Simulation time: 15.7985618114                                                   

Example 2: Use a more complicated model and more complicated delays

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.
In [3]:
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()
Info: Please mind that the delay parameters have to be set again.
Parsing file: /home/timo/Stochpy/pscmodels/CellDivision.psc
Info: No reagents have been fixed
Info: 1 trajectory is generated
simulation done!               
Info: Number of time steps 1000 End time 0.833094952035
Info: Simulation time 0.18615
Info: Delayed Next Reaction Method is selected to perform delayed stochastic simulations.
Parsing file: /home/timo/Stochpy/pscmodels/CellDivision.psc
Info: No reagents have been fixed
Info: 1 trajectory is generated
simulation done!               
Info: Number of time steps 1000 End time 0.543712164325
Info: Simulation time 0.17225
Note that the correct high-level function must be used to perform delayed stochastic simulations. We illustrate this below.
In [4]:
# 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()
Info: Delayed Direct Method is selected to perform delayed stochastic simulations.
Parsing file: /home/timo/Stochpy/pscmodels/CellDivision.psc
Info: No reagents have been fixed
*** WARNING ***: (DelayedDirectMethod) was selected. Switching to the normal Direct Method.
Info: Direct method is selected to perform stochastic simulations.
Parsing file: /home/timo/Stochpy/pscmodels/CellDivision.psc
Info: No reagents have been fixed
Info: 1 trajectory is generated
simulation done!               
Info: Number of time steps 1000 End time 0.132384906679
Info: Simulation time 0.10096
Info: Next Reaction method is selected to perform stochastic simulations
Parsing file: /home/timo/Stochpy/pscmodels/CellDivision.psc
Info: No reagents have been fixed
*** WARNING ***: an invalid method (NextReactionMethod) was selected. Switching to the Delayed Direct Method.
Info: Delayed Direct Method is selected to perform delayed stochastic simulations.
Parsing file: /home/timo/Stochpy/pscmodels/CellDivision.psc
Info: No reagents have been fixed
Info: 1 trajectory is generated
simulation done!               
Info: Number of time steps 1000 End time 1.05479142411
Info: Simulation time 0.18627
In [4]: