Simulation tools (pysb.simulator)

class pysb.simulator.CupSodaSimulator(model, tspan=None, initials=None, param_values=None, verbose=False, **kwargs)[source]

An interface for running cupSODA, a CUDA implementation of LSODA.

cupSODA is a graphics processing unit (GPU)-based implementation of the LSODA simulation algorithm (see references). It requires an NVIDIA GPU card with support for the CUDA framework version 7 or above. Further details of cupSODA and software can be found on github: https://github.com/aresio/cupSODA

The simplest way to install cupSODA is to use a pre-compiled version, which can be downloaded from here: https://github.com/aresio/cupSODA/releases

Parameters:
model : pysb.Model

Model to integrate.

tspan : vector-like, optional

Time values at which the integrations are sampled. The first and last values define the time range.

initials : list-like, optional

Initial species concentrations for all simulations. Dimensions are N_SIMS x number of species.

param_values : list-like, optional

Parameters for all simulations. Dimensions are N_SIMS x number of parameters.

verbose : bool or int, optional

Verbosity level, see pysb.simulator.base.Simulator for further details.

**kwargs: dict, optional

Extra keyword arguments, including: * gpu: Index of GPU to run on (default: 0) * vol: System volume; required if model encoded in extrinsic

(number) units (default: None)

  • obs_species_only: Only output species contained in observables (default: True)
  • cleanup: Delete all temporary files after the simulation is finished. Includes both BioNetGen and cupSODA files. Useful for debugging (default: True)
  • prefix: Prefix for the temporary directory containing cupSODA input and output files (default: model name)
  • base_dir: Directory in which temporary directory with cupSODA input and output files are placed (default: system directory determined by tempfile.mkdtemp)
  • integrator: Name of the integrator to use; see default_integrator_options (default: ‘cupsoda’)
  • integrator_options: A dictionary of keyword arguments to supply to the integrator; see default_integrator_options.

Notes

  1. If vol is defined, species amounts and rate constants are assumed to be in number units and are automatically converted to concentration units before generating the cupSODA input files. The species concentrations returned by cupSODA are converted back to number units during loading.
  2. If obs_species_only is True, only the species contained within observables are output by cupSODA. All other concentrations are set to ‘nan’.

References

  1. Nobile M.S., Cazzaniga P., Besozzi D., Mauri G., 2014. GPU-accelerated simulations of mass-action kinetics models with cupSODA, Journal of Supercomputing, 69(1), pp.17-24.
  2. Petzold, L., 1983. Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations. SIAM journal on scientific and statistical computing, 4(1), pp.136-148.
Attributes:
model : pysb.Model

Model passed to the constructor.

tspan : numpy.ndarray

Time values passed to the constructor.

initials : numpy.ndarray

Initial species concentrations for all simulations. Dimensions are number of simulations x number of species.

param_values : numpy.ndarray

Parameters for all simulations. Dimensions are number of simulations x number of parameters.

verbose: bool or int

Verbosity setting. See the base class pysb.simulator.base.Simulator for further details.

gpu : int

Index of GPU being run on

vol : float or None

System volume

n_blocks: int

Number of GPU blocks used by the simulator.

outdir : str

Directory where cupSODA output files are placed. Input files are also placed here.

opts: dict

Dictionary of options for the integrator in use.

integrator : str

Name of the integrator in use.

default_integrator_options : dict

Nested dictionary of default options for all supported integrators.
run(tspan=None, initials=None, param_values=None)[source]

Perform a set of integrations.

Returns a SimulationResult object.

Parameters:
tspan : list-like, optional

Time values at which the integrations are sampled. The first and last values define the time range.

initials : list-like, optional

Initial species concentrations for all simulations. Dimensions are number of simulation x number of species.

param_values : list-like, optional

Parameters for all simulations. Dimensions are number of simulations x number of parameters.
Returns:
A :class:`SimulationResult` object

Notes

  1. An exception is thrown if tspan is not defined in either __init__`or `run.
  2. If neither initials nor param_values are defined in either __init__ or run a single simulation is run with the initial concentrations and parameter values defined in the model.
class pysb.simulator.ScipyOdeSimulator(model, tspan=None, initials=None, param_values=None, verbose=False, **kwargs)[source]

Simulate a model using SciPy ODE integration

Uses scipy.integrate.odeint() for the lsoda integrator, scipy.integrate.ode() for all other integrators.

Warning

The interface for this class is considered experimental and may change without warning as PySB is updated.

Parameters:
model : pysb.Model

Model to simulate.

tspan : vector-like, optional

Time values over which to simulate. The first and last values define the time range. Returned trajectories are sampled at every value unless the simulation is interrupted for some reason, e.g., due to satisfaction of a logical stopping criterion (see ‘tout’ below).

initials : vector-like or dict, optional

Values to use for the initial condition of all species. Ordering is determined by the order of model.species. If not specified, initial conditions will be taken from model.initial_conditions (with initial condition parameter values taken from param_values if specified).

param_values : vector-like or dict, optional

Values to use for every parameter in the model. Ordering is determined by the order of model.parameters. If passed as a dictionary, keys must be parameter names. If not specified, parameter values will be taken directly from model.parameters.

verbose : bool or int, optional (default: False)

Sets the verbosity level of the logger. See the logging levels and constants from Python’s logging module for interpretation of integer values. False is equal to the PySB default level (currently WARNING), True is equal to DEBUG.

**kwargs : dict

Extra keyword arguments, including:

  • integrator: Choice of integrator, including vode (default), zvode, lsoda, dopri5 and dop853. See scipy.integrate.ode() for further information.
  • integrator_options: A dictionary of keyword arguments to supply to the integrator. See scipy.integrate.ode().
  • cleanup: Boolean, cleanup argument used for pysb.bng.generate_equations() call

Notes

If tspan is not defined, it may be defined in the call to the run method.

Examples

Simulate a model and display the results for an observable:

>>> from pysb.examples.robertson import model
>>> import numpy as np
>>> np.set_printoptions(precision=4)
>>> sim = ScipyOdeSimulator(model, tspan=np.linspace(0, 40, 10))
>>> simulation_result = sim.run()
>>> print(simulation_result.observables['A_total'])         
[ 1.      0.899   0.8506  0.8179  0.793   0.7728  0.7557  0.7408  0.7277
0.7158]

For further information on retrieving trajectories (species, observables, expressions over time) from the simulation_result object returned by run(), see the examples under the SimulationResult class.

run(tspan=None, initials=None, param_values=None)[source]

Run a simulation and returns the result (trajectories)

Note

tspan, initials and param_values values supplied to this method will persist to future run() calls.

Parameters:
tspan
initials
param_values

See parameter definitions in ScipyOdeSimulator.
Returns:
A :class:`SimulationResult` object
class pysb.simulator.StochKitSimulator(model, tspan=None, initials=None, param_values=None, verbose=False, **kwargs)[source]

Interface to the StochKit 2 stochastic simulation toolkit

StochKit can be installed from GitHub: https://github.com/stochss/stochkit

This class is inspired by the gillespy <https://github.com/JohnAbel/gillespy> library, but has been optimised for use with PySB.

Warning

The interface for this class is considered experimental and may change without warning as PySB is updated.

Parameters:
model : pysb.Model

Model to simulate.

tspan : vector-like, optional

Time values over which to simulate. The first and last values define the time range. Returned trajectories are sampled at every value unless the simulation is interrupted for some reason, e.g., due to satisfaction of a logical stopping criterion (see ‘tout’ below).

initials : vector-like or dict, optional

Values to use for the initial condition of all species. Ordering is determined by the order of model.species. If not specified, initial conditions will be taken from model.initial_conditions (with initial condition parameter values taken from param_values if specified).

param_values : vector-like or dict, optional

Values to use for every parameter in the model. Ordering is determined by the order of model.parameters. If passed as a dictionary, keys must be parameter names. If not specified, parameter values will be taken directly from model.parameters.

verbose : bool or int, optional (default: False)

Sets the verbosity level of the logger. See the logging levels and constants from Python’s logging module for interpretation of integer values. False is equal to the PySB default level (currently WARNING), True is equal to DEBUG.

**kwargs : dict

Extra keyword arguments, including:

  • cleanup: Boolean, delete directory after completion if True

Examples

Simulate a model and display the results for an observable:

>>> from pysb.examples.robertson import model
>>> import numpy as np
>>> np.set_printoptions(precision=4)
>>> sim = StochKitSimulator(model, tspan=np.linspace(0, 10, 5))

Here we supply a “seed” to the random number generator for deterministic results, but for most purposes it is recommended to leave this blank.

>>> simulation_result = sim.run(n_runs=2, seed=123456)

A_total trajectory for first run

>>> print(simulation_result.observables[0]['A_total'])         
[ 1.  0.  0.  0.  0.  0.]

A_total trajectory for second run

>>> print(simulation_result.observables[1]['A_total'])         
[ 1.  1.  1.  1.  0.  0.]

For further information on retrieving trajectories (species, observables, expressions over time) from the simulation_result object returned by run(), see the examples under the SimulationResult class.

run(tspan=None, initials=None, param_values=None, n_runs=1, algorithm='ssa', output_dir=None, num_processors=1, seed=None, method=None, stats=False, epsilon=None, threshold=None)[source]

Run a simulation and returns the result (trajectories)

Note

tspan, initials and param_values values supplied to this method will persist to future run() calls.

Parameters:
tspan
initials
param_values

See parameter definitions in StochKitSimulator.

n_runs : int

The number of simulation runs per parameter set. The total number of simulations is therefore n_runs * max(len(initials), len(param_values))

algorithm : str

Choice of ‘ssa’ (Gillespie’s stochastic simulation algorithm) or ‘tau_leaping’ (Tau leaping algorithm)

output_dir : str or None

Directory for StochKit output, or None for a system-specific temporary directory

num_processors : int

Number of CPU cores for StochKit to use (default: 1)

seed : int or None

A random number seed for StochKit. Set to any integer value for deterministic behavior.

method : str or None

StochKit “method” argument, default None. Only used by StochKit 2.1 (not yet released at time of writing).

stats : bool

Ask StochKit to generate simulation summary statistics if True

epsilon : float or None

Tolerance parameter for tau-leaping algorithm

threshold : int or None

Threshold parameter for tau-leaping algorithm
Returns:
A :class:`SimulationResult` object
class pysb.simulator.SimulationResult(simulator, tout, trajectories, squeeze=True, simulations_per_param_set=1)[source]

Results of a simulation with properties and methods to access them.

Warning

Please note that the interface for this class is considered experimental and may change without warning as PySB is updated.

Parameters:
simulator : Simulator

The simulator object that generated the trajectories

tout: list-like

Time points returned by the simulator (may be different from tspan if simulation is interrupted for some reason).

trajectories : list or numpy.ndarray

A set of species trajectories from a simulation. Should either be a list of 2D numpy arrays or a single 3D numpy array.

squeeze : bool, optional (default: True)

Return trajectories as a 2D array, rather than a 3d array, if only a single simulation was performed.

simulations_per_param_set : int

Number of trajectories per parameter set. Typically always 1 for deterministic simulators (e.g. ODE), but with stochastic simulators multiple trajectories per parameter/initial condition set are often desired.

Notes

In the attribute descriptions, a “trajectory set” is a 2D numpy array, species on first axis and time on second axis, with each element containing the concentration or count of the species at the specified time.

A list of trajectory sets contains a trajectory set for each simulation.

Examples

The following examples use a simple model with three observables and one expression, with a single simulation.

>>> from pysb.examples.expression_observables import model
>>> from pysb.simulator import ScipyOdeSimulator
>>> import numpy as np
>>> np.set_printoptions(precision=4)
>>> sim = ScipyOdeSimulator(model, tspan=np.linspace(0, 40, 10),                                 integrator_options={'atol': 1e-20})
>>> simulation_result = sim.run()

simulation_result is a SimulationResult object. An observable can be accessed like so:

>>> print(simulation_result.observables['Bax_c0'])         
[  1.0000e+00   1.1744e-02   1.3791e-04   1.6196e-06   1.9020e-08
   2.2337e-10   2.6232e-12   3.0806e-14   3.6178e-16   4.2492e-18]

It is also possible to retrieve the value of all observables at a particular time point, e.g. the final concentrations:

>>> print(simulation_result.observables[-1])         
(  4.2492e-18,   1.6996e-16,  1.)

Expressions are read in the same way as observables:

>>> print(simulation_result.expressions['NBD_signal'])         
[ 0.   4.7847  4.9956  4.9999  5.   5.   5.   5.   5.   5. ]

The species trajectories can be accessed as a numpy ndarray:

>>> print(simulation_result.species)
[[  1.0000e+00   0.0000e+00   0.0000e+00]
 [  1.1744e-02   5.2194e-02   9.3606e-01]
 [  1.3791e-04   1.2259e-03   9.9864e-01]
 [  1.6196e-06   2.1595e-05   9.9998e-01]
 [  1.9020e-08   3.3814e-07   1.0000e+00]
 [  2.2337e-10   4.9637e-09   1.0000e+00]
 [  2.6232e-12   6.9951e-11   1.0000e+00]
 [  3.0806e-14   9.5840e-13   1.0000e+00]
 [  3.6178e-16   1.2863e-14   1.0000e+00]
 [  4.2492e-18   1.6996e-16   1.0000e+00]]

Species, observables and expressions can be combined into a single numpy ndarray and accessed similarly. Here, the initial concentrations of all these entities are examined:

>>> print(simulation_result.all[0]) 
( 1.,  0.,  0.,  1.,  0.,  0.,  0.)

The all array can be accessed as a pandas DataFrame object, which allows for more convenient indexing and access to pandas advanced functionality, such as indexing and slicing. Here, the concentrations of the observable Bax_c0 and the expression NBD_signal are read at time points between 5 and 15 seconds:

>>> df = simulation_result.dataframe
>>> print(df.loc[5:15, ['Bax_c0', 'NBD_signal']])         
             Bax_c0  NBD_signal
time
8.888889   0.000138    4.995633
13.333333  0.000002    4.999927
all

Aggregate species, observables, and expressions trajectories into a numpy.ndarray with record-style data-type for return to the user.

dataframe

A conversion of the trajectory sets (species, observables and expressions for all simulations) into a single pandas.DataFrame.

expressions

List of trajectory sets. The first dimension contains expressions.

nsims

The number of simulations in this SimulationResult

observables

List of trajectory sets. The first dimension contains observables.

species

List of trajectory sets. The first dimension contains species.