Sphinx logo

Table Of Contents

Previous topic

StochPy’s Userguide

Next topic

The PySCeS Model Description Language

This Page

Introduction

StochPy (Stochastic Modelling in Python) is an easy-to-use package, which provides several stochastic simulation algorithms (SSA’s). These SSA’s can be used to simulate a biochemical system in a stochastic manner. Further, several unique and easy-to-use analysis techniques are provided by StochPy.

The classical approach to simulate a biochemical system is by a set of coupled ordinary differential equations (ODEs). Biological systems are often highly complex, thus there are regularly no analytical solutions available. As a result, a set of coupled ODEs is often numerically solved.

This classical approach has a deterministic nature. Therefore, a given set of ODEs will always produce the same results. These results describe the macroscopic behaviour of the biochemical system, while cell populations have regularly a heterogeneous nature. Further, the change in species amount in a deterministic model is a continuous process, while this is in reality a discrete process.

This deterministic nature becomes really problematic for systems with molecules that have a low copy number. As a result, deterministic models are often inaccurate.

SSA’s try to describe the time evolution of a reacting system, such that it takes into account discreteness and stochasticity. Discreteness and stochasticity play often important roles in biochemical systems. Therefore, SSA’s are widely used to simulate biochemical systems. Especially for biochemical systems that contain low copy numbers. For such systems, deterministic models often fail to capture the stochasticity of the system, while SSAs are capable of capturing this stochastic behaviour.

In addition, StochPy can be used as a plug-in into PySCeS. In the latest PySCeS release, a prototype is build-in to use StochPy. PySCeS - the Python Simulator for Cellular Systems - is an extendible tool kit for the analysis and investigation of cellular systems. StochPy operates completely independent from PySCeS, but it uses several parts of the PySCeS package, such as the text based Model Description Language of PySCeS. This MDL is further explained in the PySCeS Model Description Language section. PySCeS provides integrators for time simulation of deterministic models. Further, structural analysis can be done. For example, information about the stoichiometric matrix (N), the kernel matrix (K), and the link matrix (L) can be obtained. Also, metabolic control analysis (MCA) can be performed.

Start using StochPy

Here, we assume that StochPy is already installed. If not, check README.txt or installation section of this user guide.

To start modelling, it is necessary to work in an interactive Python shell. You can use iPython (recommended), IDLE (Python GUI), or simply the python command in a terminal:

$ python
Python 2.7.1+ (r271:86832, Apr 11 2011, 18:13:53)
[GCC 4.5.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>>

Then, use import stochpy:

>>> import stochpy

#######################################################################
#                                                                     #
#            Welcome to the interactive StochPy environment           #
#                                                                     #
#######################################################################
#    StochPy: Stochastic modelling in Python                          #
#    http://stompy.sourceforge.net                                    #
#    Copyright(C) T.R Maarleveld, B.G. Olivier 2010-2011              #
#    Email: tmd200@users.sourceforge.net                              #
#    VU University, Amsterdam, Netherlands                            #
#    StochPy is distributed under the BSD licence.                    #
#######################################################################

Version 1.0
Output Directory: /home/user/Stochpy
Model Directory: /home/user/Stochpy/pscmodels

This message means that the StochPy package is imported and it is ready to use. Use the next command to get more information about StochPy:

>>> help(stochpy)

One can use this help function for every available function in StochPy:

>>> mod = stochpy.SSA()
>>> help(mod)
>>> help(mod.DoStochSim)

Consequently, it is possible to start one of the four modules, which will be explained in the next sections.

Stochastic Modelling

Modelling input

The input for performing stochastic modelling are rate equations and initial conditions. For this, the PySCeS MDL is used as default format, which is further explained in the PySCeS Model Description Language section. Also, StochPy supports to use of models that are written in SBML. To use models written in SBML, several libraries are necessary, which is explained in the installation section. These SBML input files are converted into the PySCeS MDL format, which are used by the stochastic simulators.

StochPy supports assignments and SBML event facilities which resets species populations during the simulation since the StochPy 1.0.0 release.

Stochastic versus Deterministic Rate Equations

Deterministic rate equations have normally been used to describe a system of biochemical reactions, whereas these equations are often not valid for stochastic modelling. Therefore, stochastic rate equations are necessary. Therefore, it is important to understand the differences between deterministic and stochastic rate equations and why deterministic rate equations can be invalid. For this reason, the next part of this section explains the differences between deterministic and stochastic rate equations.

First, deterministic models represent species often by concentration, while stochastic models represent species often by amounts. The species amounts are identical to:

Na*[X]*Volume(L)

Here, Na is the number of Avegadro (6.02*10e23). This is the first step of the conversion of deterministic rate constants (k) to stochastic rate constants (c).

Zeroth order

Assume the following reaction:

--> X

The deterministic of rate equation of this reaction is:

k (Ms^-1)

The stochastic rate equation of this reaction is:

c (s^-1)

First order

Consider the following reaction:

X --> Y

The deterministic rate equation of this reaction is:

k*[X]

The stochastic rate equation of this reaction is:

c*x (s^-1)

Here, x corresponds to Na*[X]*Volume(L) particles. This means that x changes k*Na*[X]*Volume(L) = k*x molecules per second, which means that k=c. Thus, first order stochastic and deterministic rate equations are identical.

Second order

There exist several types of second order reactions. First, consider the next reaction:

X+Y --> Z

The deterministic rate equation of this reaction is:

k*[X]*[Y] (Ms^-1)
The stochastic rate equation of this reaction is::
c*x*y (s^-1)

Again, [X] corresponds to Na*[X]*Volume(L) particles and the same is true for [Y]. Therefore, the following relationship holds:

k*[X]*[Y] (Ms^-1) = (k*x*y)/(Na*V)^2, hence c = k/(Na*V)

Secondly, consider a dimerisation reaction:

2X --> Y

The deterministic rate equation of this reaction is:

k*X^2

The stochastic rate equation of this reaction is:

0.5*c*x*(x-1)

The first important concept to understand is that a dimerisation reaction can only occur if there are at least two molecules of species X available. For this reason, the stochastic rate equation must be zero if there is only one X molecule available.

Third order

This type of reactions does usually not occur in chemical reactions, but they are described just for the sake of completeness.

There exist three types of third order reactions:

1) X+Y+Z --> A

The deterministic rate equation of this reaction is:

k*[X]*[Y]*[Z]
The stochastic rate equation of this reaction is::

c*x*y*z

  1. X+2Y –> Z

The deterministic rate equation of this reaction is:

k*[X]*([Y]^2)

The stochastic rate equation of this reaction is:

0.5*c*x*y*(y-1)


3) 3X --> Y

The deterministic rate equation of this reaction is:

k*[X]^3

The stochastic rate equation of this reaction is:

(1/6)*x*(x-1)*(x-2)

For the last reaction, three x molecules are necessary, thus this rate equation can not fire if there are less than three molecules of x available.

Further explanations are given in the Stochastic Simulations Algorithms section.

Stochastic Simulation Algorithms

The most important module of this package is the stochastic simulation algorithm module.

Algorithms

StochPy contains four Stochastic Simulation Algorithms (SSA’s):

  • Direct Method
  • First Reaction Method
  • Next Reaction Method
  • Tau-Leaping Method

Use the following code to start using the SSA module:

>>> mod = stochpy.SSA()
Info: The Direct method is selected to perform the simulations
Parsing file: /home/user/Stochpy/pscmodels/ImmigrationDeath.psc

The default algorithm is the Direct Method and the default stochastic model is the Immigration-Death model. Next, a stochastic simulation can be done with the default settings. Of course, another method or model can be selected:

>>> mod.Method('FirstReactionMethod')
Info: The First Reaction method is selected to perform the simulations
>>> mod.Method('NextReactionMethod')
Info: The Next Reaction method is selected to perform the simulations
>>> mod.Method('TauLeaping')
Info: The Explicit Tau-Leaping method is selected to perform the simulations
>>> mod.Method('Direct')
Info: The Direct method is selected to perform the simulations

Model Selection

The following comment can be used to select another model in the PySCeS MDL:

>>> mod.Model('BurstModel.psc')
Parsing file: /home/user/Stochpy/pscmodels/BurstModel.psc
>>> mod.Model(File='BurstModel.psc',dir='/home/user/Stochpy/pscmodels')
Parsing file: /home/user/Stochpy/pscmodels/BurstModel.psc

Alternatively, one can use a model written in the systems biology markup language (SBML):

>>> mod.Model('dsmts-001-01.xml','/home/user/')
Info: extension is .xml
Info: single compartment model: locating "Death" in default compartment
Info: single compartment model: locating "Birth" in default compartment
Writing file: /home/user/Stochpy/pscmodels/dsmts-001-01.xml.psc

SBML2PSC
in : /home/user/dsmts-001-01.xml
out: /home/user/Stochpy/pscmodels/dsmts-001-01.xml.psc
Info: SBML data is converted into psc data
and is stored at: /home/user/Stochpy/pscmodels
Parsing file: /home/user/Stochpy/pscmodels/dsmts-001-01.xml.psc

>>> mod.model_file
'dsmts-001-01.xml.psc'
>>> mod.model_dir
'/home/user/Stochpy/pscmodels'

Run a SSA

Then, a simulation can be started. By default, one trajectory of 1000 time steps is generated.:

>>> mod.DoStochSim()
Info: 1 trajectory is generated
Number of time steps 1000  End time  15.2161060754
Simulation time 0.0805060863495
>>> mod.Trajectories(5)
Info: The number of trajectories is:  5
>>> mod.DoStochSim()
Info: 5 trajectories are generated
Info: Time simulation output of the trajectories is stored at ...
Info: output is written to: /home/user/Stochpy/ssa_smod.log
Simulation time 0.24516955151

Also, it is possible to select the end time of a simulation. Be careful, because this is dangerous! It is dangerous, because the reactions in the model determine when a certain reaction fires. So, dT - the time between two reactions to fire - is large if the majority of the reactions is unlikely to fire, but becomes small if the majority of the reactions is very likely to fire. As a result, model A takes 1000 time steps to reach the end time, but model B could take more than 1,000,000 time steps to reach the end time:

>>> mod.Endtime(100)

The high-level function DoStochSim() does more than starting a stochastic simulation. It accepts the following arguments:

>>> mod.DoStochSim(method = 'Direct', mode = 'time', end = 50,
trajectories = 1)
Info: The Direct method is selected to perform the simulations
Parsing file: /home/user/Stochpy/pscmodels/ImmigrationDeath.psc
Info: 1 trajectory is generated
Number of time steps 3807  End time  50.0049282696
Simulation time 0.142936944962

As a result, one can use one high-level function to determine the modelling options. Moreover, StochPy can show the the current settings:

>>> mod.ShowSpecies()
['mRNA']
>>> mod.ShowOverview()
Current Model:        ImmigrationDeath.psc
Simulation end time:  50
Current Algorithm:    <class 'stochpy.DirectMethod.DirectMethod'>
Number of trajectories:       1
Propensities are not tracked

Here, ShowSpecies() gives all the species in the model and ShowOverview() gives all the current settings.

Finally, the description of the system can be changed (the model), while you are working with StochPy. For example, you want to simulate the immigration-death model with several values of Ksyn and Kdeg. First, you start a simulation with the default values, but you decide to change these default values. Then, simply change one of the parameters in the file ImmigrationDeathmodel.psc (in the directory where the models are stored) and reload the model with the following high-level function:

>>> mod.Reload()
Parsing file: /home/user/Stochpy/pscmodels/ImmigrationDeath.psc

Build-in Analysis Techniques

After a simulation, some analysis can be done. Besides information about the time simulation, probability density functions, propensities, and waiting time plots can be created. The default line style is dotted. Furthermore, interpolation can be done if more than 1 trajectory is generated:

>>> mod.PlotTimeSim(linestyle = 'dashed',title = 'StochPy Time Simulation Plot')
>>> mod.PlotPropensities()
>>> mod.PlotDistributions()
>>> mod.PlotWaitingtimes()
>>> mod.PlotInterpolatedData()

In addition, one can use binning (default bin size = 1) for plotting of the probability density function:

>>> mod.PlotDistributions(bin_size=5)
>>> mod.PlotDistributions(bin_size=10)

Each species is plotted by default, but the user can determine which species are plotted:

>>> mod.PlotTimeSim(species2plot = ['S1','S2'])
>>> mod.PlotPropensities(rates2plot = 'R1')
>>> mod.PlotWaitingtimes(rates2plot = 'R2')
>>> mod.PlotDistributions('S4')
Error: species S4 is not in the model

In addition, stochpy.plt is available to manipulate generated plots or to make your own plots:

>>> mod.PlotTimeSim(marker = 'v')
>>> stochpy.plt.title('Your own title')
>>> stochpy.plt.xlabel('Time (s)')
>>> stochpy.plt.xlim([0,10**3])
>>> stochpy.plt.savefig('filename.pdf')   # stores the plot in the cwd

Of course, it is also possible to print such information about the simulation, which can be useful if MatPlotLib is not installed:

>>> mod.PrintTimeSim()
>>> mod.PrintPropensities()
>>> mod.PrintDistributions()
>>> mod.PrintWaitingtimes()
>>> mod.PrintMeanWaitingtimes()
>>> mod.PrintInterpolatedData()

Also, such information can be printed to a text file:

>>> mod.Write2File()
>>> mod.Write2File('TimeSim','/home/user/timesmod.txt')
>>> mod.Write2File('Propensities')
>>> mod.Write2File('Distributions')
>>> mod.Write2File('Waitingtimes')
>>> mod.Write2File('Interpol')

One can get information about the mean and standard deviation of each species during a simulation:

>>> mod.ShowMeans()
>>> mod.ShowStandardDeviations()

Finally, one can use the data objects that store the simulation data to perform their own type of analysis. See Using Stochpy as a Library.

Examples

Stochastic simulations are done with a dimerisation model to illustrate how you can use the stochastic simulation algorithm module (we assume that StochPy is already imported). This model contains two species, denoted by P and P2.

>>> mod = stochpy.SSA()
>>> mod.Model('dsmts-003-02.xml.psc')
>>> mod.TrackPropensities()
>>> mod.DoStochSim(end = 50,mode = 'time')
Info: 1 trajectory is generated
Number of time steps 591  End time  50.0230192145
Simulation time 0.0776190757751


>>> mod.PlotTimeSim()
>>> mod.PlotPropensities()
>>> mod.DoStochSim(end=5000,mode = 'time')
Info: 1 trajectory is generated
Number of time steps 17799  End time  5000.10900565
Simulation time 1.39394283295
>>> mod.PlotDistributions()
>>> mod.PlotWaitingtimes()

In addition, one can exploit the fact that StochPy is written in Python, which we illustrate with the immigration-death model. This immigration-death model contains one species, which is born born with a constant rate and dies with a first order reaction based on the amount of the species. This immigration-death model is one of the simplest examples of a Poisson process. A Poisson process is a stochastic process where events occur continuously and independently of one another.

Therefore, the mean and variance of this process are both 50. First, N samples are drawn from a Poisson distribution. Secondly, a stochastic simulation is done for N time steps. Then, the probability density function is plotted. Finally, a histogram of the randomly generated Poisson data, which can be plotted into the plot of the probability density function.:

>>> import numpy as np
>>> lambda = 50
>>> N = 2500000
>>> data = np.random.poisson(lambda,N)
>>> mod.Model('ImmigrationDeath')  # Ksyn = 10, Kdeg = 0.2, and mRNA(init) = 50
>>> mod.DoStochSim(end=N,mode='steps')
>>> mod.PlotDistributions(linestyle= 'solid')
>>> n, bins, patches = stochpy.plt.hist(data-0.5, max(data)-min(data),
    normed=1, facecolor='green')
>>> mod.ShowMeans()
mRNA  49.9117062522
>>> mod.ShowStandardDeviations()
Species       Standard Deviation
mRNA  7.0730769793

Finally, the capabilities of StochPy is illustrated with ‘burstmodel.psc’. Two different parameter sets (kon = 0.05,koff = 0.05, ksyn = 80,kdeg = 2.5 and kon = 5.0,koff = 5.0, ksyn = 80,kdeg = 2.5) are used to demonstrate this model. The mean number of mRNA will be equal, but the distribution is completely different. For the first parameter values set, the number of mRNA molecules is be bimodal distributed. In contrast, the number of mRNA molecules is unimodal distributed for the second parameter values set. Therefore, the standard deviation of the number of mRNA molecules is different. Analytical solutions are plotted in black lines.:

>>> import numpy as np
>>> import mpmath
>>> import stochpy
>>> mpmath.mp.pretty = True

>>> def GetAnalyticalPDF(kon,koff,kdeg,ksyn):
  """ Get the analytical probability density function.
  The analytical solution is taken from Sharezaei and Swain 2008
  (Analytical distributions for stochastic gene expression) """
  x_values = np.arange(0,50,0.1)
  y_values = []
  for m in x_values:
      a = ((ksyn/kdeg)**m)*np.exp(-ksyn/kdeg)/mpmath.factorial(m)
      b = mpmath.mp.gamma((kon/kdeg)+m) * mpmath.mp.gamma(kon/kdeg + koff/kdeg)/
       (mpmath.mp.gamma(kon/kdeg + koff/kdeg + m)* mpmath.mp.gamma(kon/kdeg))
      c = mpmath.mp.hyp1f1(koff/kdeg,kon/kdeg + koff/kdeg + m,ksyn/kdeg)
      y_values.append(a*b*c)
  return x_values,y_values

>>> def GetAnalyticalWaitingtimes(kon,koff,ksyn):
  """ Get analytical waiting times """
  A = mpmath.sqrt(-4*ksyn*kon+(koff + kon + ksyn)**2)
  x = []
  for i in np.arange(-20,5,0.25):
      x.append(mpmath.exp(i))
  y = []
  for t in x:
      B = koff + ksyn - (mpmath.exp(t*A)*(koff+ksyn-kon))-kon+A+ mpmath.exp(t*A)*A
      p01diff = mpmath.exp(-0.5*t*(koff + kon + ksyn+A))*B/(2.0*A)
      y.append(p01diff*ksyn)
  return (x,y)

  >>> mod = stochpy.SSA()
  >>> mod.Model('Burstmodel.psc')  # Parameter values in Burstmodel.psc: kon = koff = 0.05
  >>> ntimesteps = 1000000
  >>> mod.DoStochSim(end=ntimesteps,mode='steps')
  >>> mod.PlotDistributions(species2plot = 'mRNA', colors=['#32CD32'],marker = 'o') # Usage of html color codes
  >>> mod.PlotWaitingtimes('R3', colors=['#32CD32'],linestyle = 'None',marker='o')

  >>> ### Change Parameter values in Burstmodel.psc: kon = koff = 5.0
  >>> mod.Reload()
  >>> mod.DoStochSim(end=ntimesteps,mode='steps')
  >>> stochpy.plt.figure(1)
  >>> mod.PlotDistributions(species2plot = 'mRNA', colors=['r'],marker ='v')

  >>> kon = 0.05
  >>> koff = 0.05
  >>> kdeg = 2.5
  >>> ksyn = 80.0
  >>> x,y = GetAnalyticalPDF(kon,koff,kdeg,ksyn)
  >>> stochpy.plt.step(x,y,color ='k')

  >>> kon = 5.0
  >>> koff = 5.0
  >>> x,y = GetAnalyticalPDF(kon,koff,kdeg,ksyn)
  >>> stochpy.plt.step(x,y,color ='k')
  >>> stochpy.plt.xlabel('mRNA copy number per cell')
  >>> stochpy.plt.ylabel('Probability mass')
  >>> stochpy.plt.legend(['Bimodal','Unimodal', 'Analytical solution'],numpoints=1,frameon=False)
  >>> stochpy.plt.title('')
  >>> stochpy.plt.ylim([0,0.045])

  >>> ##################### Finish the waiting times plot #######################

  >>> ### Change Parameter values ##
  >>> mod.PlotWaitingtimes('R3', colors=['r'],linestyle = 'None',marker='v')

  >>> kon = 0.05
  >>> koff = 0.05
  >>> (x,y) = GetAnalyticalWaitingtimes(kon,koff,ksyn)
  >>> stochpy.plt.plot(x,y,color ='k')

  >>> kon = 5.0
  >>> koff = 5.0
  >>> (x,y) = GetAnalyticalWaitingtimes(kon,koff,ksyn)
  >>> stochpy.plt.plot(x,y,color ='k')
  >>> stochpy.plt.xlabel('Time between RNA synthesis events')
  >>> stochpy.plt.ylabel('Probability density')
  >>> stochpy.plt.legend(['Bimodal','Unimodal', 'Analytical solution'],numpoints=1,frameon=False,loc='lower left')
  >>> stochpy.plt.title('')
  >>> stochpy.plt.xlim([10**-7,10**3])
  >>> stochpy.plt.ylim([10**-9,10**3])

Using StochPy as a Library

It is straightforward to use StochPy as a library in your code. A data object, data_stochsim, is created for those that want to use StochPy as a library or for people that want to do their own analysis. Species, distributions, propensities, simulation time, and waiting times are stored in for instance NumPy arrays and lists. In addition, labels are stored for each of these data types in separate lists. Of course, data such as distributions are only available if they are calculated.

Furthermore, determined means and standard deviations are stored in dictionaries. Finally, information about the stochastic simulation such as the number of time steps, the simulation end time, and the trajectory is stored.

This data object (data_stochsim) is written to disk space if multiple trajectories are generated. The high-level function GetTrajectoryData(n) can be used to get access to the simulation data of a specific trajectory. By default, the latest generated trajectory is not written to disk space, thus accessible without using the high-level function GetTrajectoryData(n):

>>> import stochpy
>>> mod = stochpy.SSA()
>>> mod.DoStochSim(trajectories = 10,mode = 'steps',end = 1000)
>>> mod.data_stochsim
<stochpy.PyscesMiniModel.IntegrationStochasticDataObj object at 0x32cd750>
>>> mod.data_stochsim.simulation_trajectory
10
>>> mod.data_stochsim.time          # time array (not shown)
>>> mod.data_stochsim.species       # species array (not shown)
>>> mod.data_stochsim.species_labels
>>> mod.data_stochsim.getSpecies()  # time + species array (not shown)
>>> mod.GetMeans()                  # for each species
>>> mod.data_stochsim.means
>>> mod.GetDistributions()
>>> mod.data_stochsim.distributions # for each species (not shown)
>>> mod.GetWaitingtimes()
>>> mod.data_stochsim.waiting_times # for each reaction (not shown)
>>> mod.GetTrajectoryData(5)
>>> mod.data_stochsim.simulation_trajectory
5

Alternatively, one can be interested in interpolated data from one simulation with multiple trajectories. A second data object (data_stochsim_interpolated) is becomes available if the user uses one of the functions to get interpolated data:

>>> mod.data_stochsim_interpolated
AttributeError: SSA instance has no attribute 'data_stochsim_interpolated'
>>> mod.GetInterpolatedData()
>>> mod.data_stochsim_interpolated
<stochpy.PyscesMiniModel.InterpolatedDataObj object at 0x32cd3d0>
>>> mod.data_stochsim_interpolated.time  # time array (not shown)
>>> mod.data_stochsim_interpolated.means # at every t (not shown)
>>> mod.data_stochsim_interpolated.standard_deviations  # SDs (not shown)

Stochastic Test suite

The stochastic test suite from Evans et al. 2008 (The SBML discrete stochastic models test suite) is used to test StochPy. This test suite tests stochastic simulation software on the following points:

  • local and global parameters (parameter overloading)
  • boundary conditions
  • Cell compartment volume
  • hasOnlySubstanceUnits flag
  • math expression parsing
  • compartment volume explicitly including in the rate laws
  • assignment rules
  • time events
  • species population events

StochPy successfully reproduces the desired results, whereas there is one exception. StochPy converts species concentrations (HasOnlySubstanceUnits = True) to species amounts. Here, the species concentration is multiplied by the volume of the compartment. In the stochastic test suite results, a simulation is done with species concentrations, which gives a different simulation result.

Some examples (dsmts-001-07, dsmts-001-19, dsmts-003-03, and dsmts-003-04) are shown here. These models test on multiple species, assignment rules, time events, and species amount events respectively.

Nucleosome Modification Simulations

As mentioned in the Stochastic Simulation Algorithms section StochPy contains four stochastic simulation algorithms. These algorithms can be used to perform simulations on all sorts of biochemical systems.

Usually, default time plots, distributions, propensities, and waiting times are sufficient analysis techniques. However, this is not always to case. Here, we demonstrate the flexibility of StochPy, because it is relatively easy to add your own modules or analysis techniques.

One example is the simulation of nucleosome modification models. Here, each nucleosome can carry several modifications. As a result, each nucleosome is described through several reactions, where each reaction gives information about a particular modification. Here, we are interested in for example the global pattern of these modifications. Therefore, default analysis techniques are insufficient. For this reasons, a nucleosome simulation module was built, which can be used interactively (make sure that StochPy is imported, which is shown in the Start using StochPy section ):

>>> smod = stochpy.NucSim()
Welcome to the nucleosome modification simulation module
Info: The Direct method is selected to perform the simulations
Parsing file: /home/user/Stochpy/pscmodels/model1.psc

Now, it is possible to start a simulation with the default settings, but it is again possible to select for example another model or choose another number of time steps, just as is described in the Stochastic Simulation Algorithms section:

>>> smod.GetGapMeasure()
Gap Measure   0.993
>>> smod.PlotGlobalTimeSim()
>>> smod.PlotGlobalDistributions()
>>> smod.PlotPattern()
>>> smod.PlotStateTimes('M1')

Again, it is possible to plot information from not all species:

>>> smod.PlotGlobalTimeSim('M')
>>> smod.PlotGlobalTimeSim(['M','A'])

And to print the output or write it to a file:

>>> smod.PrintGlobalTimeSim()
>>> smod.PrintGlobalDistributions()
>>> smod.PrintPattern()
>>> smod.PrintStateTimes()
>>> smod.Write2File()

Nucleosome Model Builder

An average gene contains about 40 nucleosomes. Therefore, nucleosome modification models get enormously large. As a result, a nucleosome model builder module is created, which quickly builds the user-defined models:

>>> model = stochpy.NucModel()
>>> model.Build()
/home/user/Stochpy/pscmodels
The generated model is stored at /home/user/Stochpy/pscmodels/model1.psc
>>> help(model)

By default, a nucleosome model with 20 nucleosomes was built which can be considered as a small gene. A small gene was chosen, because the number of reactions explodes for larger number of nucleosomes and 20 nucleosomes was still large enough to get all sorts of biological phenomena. Further, 8 pre-defined models are build-in:

>>> model.ModelType(2)
ModelType:    2

Of course, the number of nucleosomes can be changed:

>>> model.ChangeN(10)
The number of nucleosomes is:         10
>>> model.Build()
/home/user/Stochpy/pscmodels
The generated model is stored at /home/user/Stochpy/pscmodels/model2.psc

Notice that the pre-defined models are build at the configuration step of StochPy, thus these are overwritten if new versions is build.

Also, landing zones interactively can be determined interactively:

>>> model.BuildLandingZones({'M',[10,11]})

Finally, several other high-level functions are developed:

>>> model.Recruitment()
Recruitment is activated
>>> model.Neighbours()
Info: Neighbour interactions are activated
>>> model.LongRangeInteractions()
Info: Long range neighbour interactions are activated"
Info: The default threshold is 7.
      Use Threshold(value) to change the threshold.
>>> model.ChangeThreshold(5)
Info: The Threshold is        5
>>> model.Decay()
Info: (Long range) neighbour interactions are activated
Info: The default threshold is 7.
      Use Threshold(value) to change the threshold.
Info: Decay effect is activated

Installation and Configuration

Installation

The following software is required before installing StochPy:

The Following software is optional but recommended:

  • Matplotlib
  • libsbml
  • libxml2

libsbml is necessary to convert models written in SBML format to the PySCeS MDL. The libsbml software requires a XML parser library which is libxml2 by default.

Windows

Use the windows installer (StochPy-1.0.0-linux-i686.exe) in directory: /StochPy-1.0.0/

Linux/MAC OS/Cygwin

In the command line, as root (sudo -s):

$ cd /.../StochPy-1.0.0
$ python setup install

Configuration

The StochPy package contains some example models, such as:

  • Immigration-Death model
  • Burst model
  • Decaying-Dimerizing model
  • Prokaryotic auto-regulation model

All these files are placed in the directory /home dir/Stochpy/pscmodels/ after the first usage of the package. All these models are written in the PySCeS MDL. This PySCeS MDL is explained in the PySCeS Model Description Language section. Of course, SBML models can be used as input. A SBML2PySCeS format converter is available which is used automatically if a SBML file is given as input.