Examples

Click here to download the StochPy TestScript to do some basic examples. Alternatively, one can exploit the fact that StochPy is written in Python, which we illustrate with the following examples. All these examples are hard coded in the StochPy Utilities module, which is available since the release of 1.0.2. In these examples, most of the plots start as a default StochPy plot, whereas they are altered to the requirements of the user. This can be done with basic commands as is illustrated.

In future releases, more examples and useful functions built by users of StochPy will be added to the utils module.

Usage of this module is done with the following code:

import stochpy
utils = stochpy.Utils()
utils.doExample1()

1) Immigration-Death Model (ImmigrationDeath.psc)

The immigration-death model contains one species, which is born with a constant rate and dies with a first order reaction based on the copy number of that species. Because StochPy is designed to model biological systems, we assume constant syntheses of mRNA molecules and a first order degradation of mRNA molecules. 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 copy numbers per cell. 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 is plotted on top of the probability density function of the mRNA copy number in the cell. A cell size of 1 was taken for convenience.

>>> import stochpy
>>> mod = stochpy.SSA()
>>> import numpy as np
>>> lambda = 50
>>> N = 1500000
>>> 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

2) Dimerisation Model (dsmts-003-04.psc)

The dimerisation model contains two species, denoted by P and P2, representing a dimerisation process. The initial numbers of molecules of P and P2 are 100 and 0, respectively. A cell size of 1 was taken for convenience.
The deterministic rate equation of this reaction is: k*P^2, while the stochastic rate equation is 0.5*P*(P-1). This dimerisation reaction can only occur if there are at least two molecules of species P available in the cell. For this reason, the stochastic rate equation must be zero if there is only one P molecule available in the cell.
In this specific example, the species populations in the cell are reset to P = 100 and P2 = 0 if the number of P2 molecules becomes greater than 30 in the cell. This example is taken from the Stochastic Test Suite of Evans et al. 2008.

>>> mod = stochpy.SSA()
>>> mod.Model('dsmts-003-04.psc')
>>> mod.DoStochSim(end = 50,mode = 'time',trajectories = 10000)
>>> mod.PlotInterpolatedData()

images/dsmts-003-04-PlotInterpolatedData.svg

3) Burst Model (Burstmodel.psc)

The model shown in the figure above describes stochastic single-cell transcription. This transcription can occur in a bursty and non-bursty manner, which depends on the used parameter values. Two different parameter sets (kon = 0.05 per min, koff = 0.05 per min, ksyn = 80 per min, kdeg = 2.5 per min and kon = 5.0 per min, koff = 5.0 per min, ksyn = 80 per min, kdeg = 2.5 per min) are used to demonstrate this model. A cell volume size of 1 was taken for convenience.

Part I

First, two "short" time simulations were done to illustrate the behaviour of mRNA through time for bursty and non-bursty transcription. HTML color codes were used to plot species with specific colors.

>>> import stochpy
>>> import matplotlib.gridspec as gridspec
>>> sim_end = 100
>>> mod = stochpy.SSA()
>>> mod.Model('Burstmodel.psc')
>>> mod.DoStochSim(end=sim_end,mode='time',trajectories = 1)
>>> # Use a nice grid to plot 4 figures
>>> gs = gridspec.GridSpec(4,1,width_ratios=[1],height_ratios=[0.3,1,0.3,1])
>>> ax1 = stochpy.plt.subplot(gs[0])
  
>>> mod.PlotTimeSim('ONstate')
>>> stochpy.plt.xlabel('')               # remove xlabel
>>> stochpy.plt.ylabel('')               # remove ylabel
>>> stochpy.plt.legend('',frameon=False) # remove legend
>>> stochpy.plt.xlim([0,sim_end])        # set x lim
>>> stochpy.plt.xticks([])               # remove x ticks
>>> stochpy.plt.ylim([0,1.1])            # set y lim
>>> stochpy.plt.yticks([])               # remove y lim
>>> stochpy.plt.text(-5.5,0.9,'ON')      
>>> stochpy.plt.text(-5.5,0,'OFF')
>>> stochpy.plt.text(101,0.35,'A',fontsize = 14)

>>> ax2 = stochpy.plt.subplot(gs[1])
>>> mod.plot.ResetPlotnum()
>>> mod.PlotTimeSim('mRNA',colors = ['#32CD32'])
>>> stochpy.plt.xlim([0,sim_end])
>>> stochpy.plt.legend('',frameon=False)
>>> stochpy.plt.xticks([])
>>> stochpy.plt.title('')
>>> stochpy.plt.xlabel('')
>>> stochpy.plt.ylabel('mRNA')
>>> stochpy.plt.yticks([0,20,40,60])
>>> stochpy.plt.text(101,27,'B',fontsize = 14)

>>> a = raw_input('Change parameters kon and koff to 5.0 and press any key to continue')

>>> mod.Reload()

>>> ax3 = stochpy.plt.subplot(gs[2])
>>> mod.plot.ResetPlotnum()
>>> mod.DoStochSim(end=sim_end,mode='time')
>>> mod.PlotTimeSim('ONstate')
>>> stochpy.plt.title('')
>>> stochpy.plt.xlabel('')
>>> stochpy.plt.ylabel('')
>>> stochpy.plt.legend('',frameon=False)
>>> stochpy.plt.xlim([0,sim_end])
>>> stochpy.plt.xticks([])
>>> stochpy.plt.ylim([0,1.1])
>>> stochpy.plt.yticks([])
>>> stochpy.plt.text(-5.5,0.9,'ON')
>>> stochpy.plt.text(-5.5,0,'OFF')
>>> stochpy.plt.text(101,0.35,'C',fontsize = 14)

>>> ax4 = stochpy.plt.subplot(gs[3])
>>> mod.plot.ResetPlotnum()
>>> mod.PlotTimeSim('mRNA',colors = ['r'])
>>> stochpy.plt.xlim([0,sim_end])
>>> stochpy.plt.legend('',frameon=False)
>>> stochpy.plt.xticks([])
>>> stochpy.plt.title('')
>>> stochpy.plt.ylabel('mRNA')
>>> stochpy.plt.yticks([0,20,40,60])
>>> stochpy.plt.text(101,27,'D',fontsize = 14)
>>> stochpy.plt.savefig('stochpy_test.pdf')


 Part II 
         
         The mean copy number of mRNA in the cell is identical for bursty and non-bursty transcription, but the mRNA copy number distribution is completely different. For the first parameter values set, the mRNA copy number is bimodal distributed. In contrast, the mRNA copy number is unimodal distributed for the second parameter values set. Therefore, the standard deviation of the mRNA copy number is different. 
>>> 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 = 1500000
>>> mod.DoStochSim(end=ntimesteps,mode='steps')
>>> mod.ShowMeans()
mRNA 16.2853951455 
ONstate 0.508926320062 
OFFstate 0.491073679938

>>> mod.ShowStandardDeviations()
Species Standard Deviation
mRNA  16.1898452499 
ONstate 0.49992031446 
OFFstate 0.49992031446 
>>> mod.PlotDistributions(species2plot = 'mRNA', colors=['#00FF00'],linestyle = 'solid')
>>> mod.PlotWaitingtimes('R3', colors=['##00FF00'],linestyle = 'solid')
 
                 Change parameter values in Burstmodel.psc: kon = koff = 5.0
 
>>> mod.Reload()
>>> mod.DoStochSim(end=ntimesteps,mode='steps')
>>> mod.plot.plotnum=1
>>> mod.PlotDistributions(species2plot = 'mRNA', colors='k',linestyle = 'solid')
>>> mod.ShowMeans()
mRNA 16.018415322 
ONstate 0.500097112814 
OFFstate 0.499902887186
>>> mod.ShowStandardDeviations()
Species       Standard Deviation
mRNA  8.2433278546 
ONstate 0.499999990569
OFFstate 0.499999990569
 
                 Get the analytical probability distribution for two different sets of parameter values. This analytical solution is taken from Sharezaei and Swain 2008 (Analytical distributions for stochastic gene expression). 

>>> kon = 0.05
>>> koff = 0.05
>>> 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.legend(['Bimodal','Unimodal', 'Analytical solution'],numpoints=1,frameon=False)

                 
                 Finish the waiting times plot by plotting the by StochPy calculated waiting times and two analytical solution for two different 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.figure(2)
>>> >>> 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])