StochPy SSA Module Example

We use this iPython notebook to illustrate the basics of StochPy. Written by TR Maarleveld, Amsterdam, The Netherlands E-mail: [email protected] Last Change: August 10, 2015

In [2]:
import stochpy
smod = stochpy.SSA(IsQuiet=True)
# required for iPython Notebook inline plotting
%matplotlib inline 
#######################################################################
#                                                                     #
#            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-2015 #
#  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.3.0
Output Directory: /home/timo/Stochpy
Model Directory: /home/timo/Stochpy/pscmodels

1: Basic Simulation with the Direct method

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.
In [3]:
smod.DoStochSim(IsTrackPropensities=True)
smod.PlotSpeciesTimeSeries()
smod.PlotPropensitiesTimeSeries()
In [6]:
smod.GetWaitingtimes()
smod.PlotWaitingtimesDistributions()
smod.PrintWaitingtimesMeans()
stochpy.plt.xlim(xmin=10**-3)
Reaction	Mean
R1	0.098
R2	0.101
Out[6]:
(0.001, 1.0)

2: Export data to a text file

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).
In [6]:
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')
*** WARNING ***: Autocorrelations are not yet calculated. StochPy automatically calculates autocorrelations with pre-defined settings. You can use GetSpeciesAutocorrelations(species2calc=True,n_samples=51)

3: More stochastic simulations

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
In [7]:
# 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()
3
Species	Mean
mRNA	50.286
Species	Standard Deviation
mRNA	6.773
In [8]:
# Switch to data from trajectory 1 and show the means of each species
smod.GetTrajectoryData(1)
smod.PrintSpeciesMeans()
smod.PrintSpeciesStandardDeviations()
Species	Mean
mRNA	50.205
Species	Standard Deviation
mRNA	6.109
In [9]:
# Do one long simulation
smod.DoStochSim(trajectories=1,end=10**6,mode='steps')
smod.PrintSpeciesMeans()
smod.PrintSpeciesStandardDeviations()
Species	Mean
mRNA	50.086
Species	Standard Deviation
mRNA	7.123
In [11]:
# Plot the PDF for different bin sizes
smod.PlotSpeciesDistributions()
smod.PlotSpeciesDistributions(bin_size=5)  # larger bin size
smod.Export2File(analysis='distribution',datatype='species')

4: Algorithms

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
In [13]:
smod.Model('DecayingDimerizing.psc')
smod.DoStochSim(method='tauleaping',end=50,mode='time')
smod.PlotSpeciesTimeSeries()
stochpy.plt.xscale('log')

5: Model Format

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.
In [15]:
# 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
Species	Mean
mRNA	99.434
In [16]:
# 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)
Info: single compartment model: locating "Death" in default compartment
Info: single compartment model: locating "Birth" in default compartment
Writing file: /home/timo/Stochpy/pscmodels/dsmts-001-01.xml.psc
                                                   
In [20]:
# 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()
csymbol time defined as "t" in event 
Info: single compartment model: locating "Dimerisation" in default compartment
Info: single compartment model: locating "Disassociation" in default compartment
Writing file: /home/timo/Stochpy/pscmodels/dsmts-003-03.xml.psc
                                                   
In [21]:
# 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()
                                                   
In [22]:
# Volume Models
smod.Model('dsmts-001-11.xml.psc') 
smod.DoStochSim(method = 'Direct',trajectories=100,end=50,mode ='time')
smod.GetRegularGrid()
smod.PlotAverageSpeciesTimeSeries()
*** WARNING ***: Species are given in concentrations and the compartment volume is unequal to one.
The species concentrations are multiplied by the compartment volume.
*** WARNING ***: Species are given in concentrations and the compartment volume is unequal to one.
The species concentrations are multiplied by the compartment volume.
                                                   
In [ ]: