run_ssa in PySB library working example?

Question

I tried run_ssa from Python PySB library using "Hello world" ODE simulation model (replacing the concentrations with corresponding molecular number and converting deterministic rates to stochastic equivalents).

This is trajectories obtained using odesolve method:

[(301.0, 180.0, 0.0, 0.0)
 (269.7114893451928, 148.71148934519294, 31.288510654807023, 31.288510654807023)
 (246.1336243593057, 125.13362435930587, 54.866375640694116, 54.866375640694116)
 (227.7848590843455, 106.78485908434567, 73.21514091565433, 73.21514091565433)
 (213.1444366774497, 92.1444366774498, 87.8555633225502, 87.8555633225502)...

which look reasonable.

And this is what run_ssa returns:

[(0.0, 0.0) (0.4, 27.0) (0.8, 54.0) (1.2, 73.0) (1.6, 81.0) (2.0, 98.0)
 (2.4, 116.0) (2.8, 124.0) (3.2, 129.0) (3.6, 135.0) (4.0, 137.0)
 (4.4, 141.0) (4.8, 144.0) (5.2, 145.0) (5.6, 146.0) (6.0, 149.0)
 (6.4, 153.0) (6.8, 157.0) (7.2, 159.0) (7.6, 160.0) (8.0, 162.0)
 (8.4, 162.0) (8.8, 165.0) (9.2, 164.0) (9.6, 166.0) (10.0, 167.0)
 (10.4, 168.0) (10.8, 171.0) (11.2, 172.0) (11.6, 172.0) (12.0, 174.0)
 (12.4, 175.0) (12.8, 176.0) (13.2, 176.0) (13.6, 176.0) (14.0, 176.0)
 (14.4, 177.0) (14.8, 177.0) (15.2, 177.0) (15.6, 177.0) (16.0, 177.0)
 (16.4, 177.0) (16.8, 177.0) (17.2, 177.0) (17.6, 177.0) (18.0, 177.0)
 (18.4, 177.0) (18.8, 177.0) (19.2, 177.0) (19.6, 178.0) (20.0, 178.0)
 (20.4, 178.0) (20.8, 178.0) (21.2, 178.0) (21.6, 178.0) (22.0, 178.0)
 (22.4, 178.0) (22.8, 179.0) (23.2, 179.0) (23.6, 179.0) (24.0, 180.0)
 (24.4, 180.0) (24.8, 180.0) (25.2, 180.0) (25.6, 180.0) (26.0, 180.0)
 (26.4, 180.0) (26.8, 180.0) (27.2, 180.0) (27.6, 180.0) (28.0, 180.0)
 (28.4, 180.0) (28.8, 180.0) (29.2, 180.0) (29.6, 180.0) (30.0, 179.0)
 (30.4, 179.0) (30.8, 179.0) (31.2, 179.0) (31.6, 179.0) (32.0, 179.0)
 (32.4, 179.0) (32.8, 179.0) (33.2, 179.0) (33.6, 179.0) (34.0, 179.0)
 (34.4, 179.0) (34.8, 179.0) (35.2, 180.0) (35.6, 180.0) (36.0, 180.0)
 (36.4, 180.0) (36.8, 180.0) (37.2, 180.0) (37.6, 180.0) (38.0, 180.0)
 (38.4, 180.0) (38.8, 180.0) (39.2, 180.0) (39.6, 180.0) (40.0, 180.0)]

Looks like instead of taking the initial conditions as a start point, it uses one of them as a limiting factor? And the first array looks like just time to me... According to description it is supposed to return trajectories just as ODE simulation...

By the look of it odesolve calles python library to solve ODE but run_ssa calls BioNetGen modelling language.

I've tried to implement both deterministic and stochastic models for Michaelis-Menten kinetics directly using BioNetGen and they looked fine.

So, the question is: has anyone ever got it working? Any working examples of using it?

The purpose of this is testing the library and I just need to know if it is possible to make it work and if anyone ever managed to do so.

Answer 1

PySB author/developer here.

It appears that the fault lies with us for (I just discovered after reading this SO post) not documenting the return type of the run_ssa function. This will be remedied. As mentioned by AliceNCL, the output returned by run_ssa is parsed back from the BNG .gdat file. However, the trajectory returned is not random: the .gdat file contains the trajectory for the only Observable that was defined in the model: 'LR', the ligand-receptor complex. The .cdat file, on the other hand, contains trajectories for all of the species defined in the model, simply numbered 1, 2, 3. So, as currently implemented, if there is a particular species you wish to get the SSA trajectory for, you would add an observable for it.

The odesolve function on the other hand, returns a recarray with all of the species (indexed as '_s0', '_s1', '__s2') AND the named observable 'LR'. Hence, there are four entries, and the last two are the same: the ligand-receptor complex as "species 2" and as the named observable "LR".

If you (or anyone else) has any other issues using the library, please let us know by emailing the PySB mailing list (see link at pysb.org) or posting an issue on GitHub! We check Stack Overflow only very sporadically.

In any case, here is a working example:

from pysb import *
from pysb.examples.hello_pysb import model
from pylab import *
from pysb import bng
from pysb.integrate import odesolve

ion()
t_end = 100

# BNG SSA simulation
ssa_sim = bng.run_ssa(model, t_end=t_end, n_steps=1000)
figure()
plot(ssa_sim['time'], ssa_sim['LR'], color='r')

# ODE simulation
t = linspace(0, t_end)
ode_sim = odesolve(model, t)
plot(t, ode_sim['LR'], color='g')

which produces the following: https://i.stack.imgur.com/wPtOH.jpg

Answer 2

Thanks, I did realise what exactly it actually does in the end and managed to produce some working examples. But it took some digging through the code and BNG documentation :). The thing is that odesolve in the documentation example actually returns the whole lot of trajectories and you choose the ones you want when plotting, not just observable ones so it got me confused. Certainly, the addition of this working example and some more explanation to the documentation will be definitely a great help to the researches! Other then this it is a good and useful library