We use this iPython notebook to illustrate the basics of StochPy.

Written by TR Maarleveld, Amsterdam, The Netherlands
E-mail: tmd200@users.sourceforge.net
Last Change: August 10, 2015 
import stochpy
smod = stochpy.SSA(IsQuiet=True)
# required for iPython Notebook inline plotting
%matplotlib inline 
We start by using the direct method for basic simulation of the immigration-death model (the default model). The immigration-death model is a simplification of the birth-death model (a Markovian process) with a fixed birth rate.

We generate a complete and detailed history of all species in the model (mRNA), the propensities of the immigration and death reaction, and which reaction fires when. This allows us to determine the event waiting times for each event, i.e. reaction.

This is what we illustrate here.
smod.DoStochSim(IsTrackPropensities=True)
smod.PlotSpeciesTimeSeries()
smod.PlotPropensitiesTimeSeries()

smod.GetWaitingtimes()
smod.PlotWaitingtimesDistributions()
smod.PrintWaitingtimesMeans()
stochpy.plt.xlim(xmin=10**-3)
StochPy also allows the user to export the data to a text file via the high-level function "Export2File".
By default, the species time series data is exported, whereas all types of data can be exported (if generated of course).
smod.Export2File()
smod.Export2File(analysis='distribution')
smod.Export2File(analysis='distribution',datatype='species')
smod.Export2File(analysis='mean',datatype='species')
smod.Export2File(analysis='std',datatype='species')
smod.Export2File(analysis='autocorrelation',datatype='species')
We again use to immigration-death model to:
- generate multiple realizations (trajectories)
- show how we can switch between the different generated trajectories within StochPy
- perform a simulation for 10**6 time steps to illustrate probability density functions of species copy numbers
# Show the means from the data of last generated trajectory
smod.Model('ImmigrationDeath.psc')
smod.DoStochSim(trajectories=3,end=10**3) # multiple trajectories
print(smod.data_stochsim.simulation_trajectory)
smod.PrintSpeciesMeans()
smod.PrintSpeciesStandardDeviations()

# Switch to data from trajectory 1 and show the means of each species
smod.GetTrajectoryData(1)
smod.PrintSpeciesMeans()
smod.PrintSpeciesStandardDeviations()


# Do one long simulation
smod.DoStochSim(trajectories=1,end=10**6,mode='steps')
smod.PrintSpeciesMeans()
smod.PrintSpeciesStandardDeviations()

# Plot the PDF for different bin sizes
smod.PlotSpeciesDistributions()
smod.PlotSpeciesDistributions(bin_size=5)  # larger bin size
smod.Export2File(analysis='distribution',datatype='species')
StochPy provides many different stochastic simulation algorithms which we can divide in the following categories:
1) exact SSAs
2) inexact SSA
3) exact SSAs with delay
4) single molecule methods

The exact SSAs are most used, but suffer typically from extensive running times. Inexaxt SSAs (i.e. tau-leaping) can be an outcome when exact SSAs are too slow for the job. We illustrate this here for the "DecayingDimerizing" model. The simulation with an exact solver is about 50 times slower. Remember that using tau-leaping methods result in loss of information (e.g. waiting times cannot be calculated).

Exact SSAs with delay and single molecule methods are discussed in a different iPython notebook
smod.Model('DecayingDimerizing.psc')
smod.DoStochSim(method='tauleaping',end=50,mode='time')
smod.PlotSpeciesTimeSeries()
stochpy.plt.xscale('log')
StochPy supports e.g.
- models written in the PySCeS MDL (human readable) and SBML. 
- assignments and events (time and species)
- modification of model parameters and species copy number amounts

In the upcoming examples, we illustrate each of these features in more detail.
# Model modification
smod = stochpy.SSA() # loads the default immigration-death model
smod.ChangeParameter('Ksyn',20.0)
smod.ChangeParameter('Kdeg',0.2)
smod.ChangeInitialSpeciesCopyNumber('mRNA',100)
smod.DoStochSim(end=10**5)
smod.PrintSpeciesMeans()   # should be ~Ksyn/Kdeg

# Import a model written in SBML and generate a regular grid to compare different trajectories
smod.Model('dsmts-001-01.xml')
smod.DoStochSim(trajectories=100,end=50,mode='time') 
smod.GetRegularGrid()
smod.PlotAverageSpeciesTimeSeries()
smod.Export2File(datatype='species',analysis='timeseries',IsAverage=True)

# Test a model with a time event
smod.Model('dsmts-003-03.xml.psc') 
smod.DoStochSim(trajectories=1000,end=50,mode='time')
smod.GetRegularGrid()
smod.PlotAverageSpeciesTimeSeries()

# Use the First Reaction method to test a model with a species amount event 
smod.Model('dsmts-003-04.xml.psc') 
smod.DoStochSim(method = 'FirstReactionMethod',trajectories=1000,end=50,mode='time')
smod.GetRegularGrid()
smod.PlotAverageSpeciesTimeSeries()

# Volume Models
smod.Model('dsmts-001-11.xml.psc') 
smod.DoStochSim(method = 'Direct',trajectories=100,end=50,mode ='time')
smod.GetRegularGrid()
smod.PlotAverageSpeciesTimeSeries()