Simulation tools (`pysb.simulator`)¶

class pysb.simulator.BngSimulator(model, tspan=None, cleanup=True, verbose=False)[source]¶

Simulate a model using BioNetGen

run(tspan=None, initials=None, param_values=None, n_runs=1, method='ssa', output_dir=None, output_file_basename=None, cleanup=True, population_maps=None, **additional_args)[source]¶

Simulate a model using BioNetGen

Parameters:

tspan: vector-like

time span of simulation

initials: vector-like, optional

initial conditions of model

param_values : vector-like or dictionary, optional

Values to use for every parameter in the model. Ordering is determined by the order of model.parameters. If not specified, parameter values will be taken directly from model.parameters.

n_runs: int

number of simulations to run

method : str

Type of simulation to run. Must be one of:

‘ssa’ - Stochastic Simulation Algorithm (direct method with propensity sorting)

‘nf’ - Stochastic network free simulation with NFsim. Performs Hybrid Particle/Population simulation if population_maps argument is supplied

‘pla’ - Partioned-leaping algorithm (variant of tau-leaping algorithm)

‘ode’ - ODE simulation (Sundials CVODE algorithm)

output_dir : string, optional

Location for temporary files generated by BNG. If None (the default), uses a temporary directory provided by the system. A temporary directory with a random name is created within the supplied location.

output_file_basename : string, optional

This argument is used as a prefix for the temporary BNG output directory, rather than the individual files.

cleanup : bool, optional

If True (default), delete the temporary files after the simulation is finished. If False, leave them in place. Useful for debugging.

population_maps: list of PopulationMap

List of PopulationMap objects for hybrid particle/population modeling. Only used when method=’nf’.

additional_args: kwargs, optional

Additional arguments to pass to BioNetGen

Examples

Simulate a model using network free simulation (NFsim):

>>> from pysb.examples import robertson
>>> from pysb.simulator.bng import BngSimulator
>>> model = robertson.model
>>> sim = BngSimulator(model, tspan=np.linspace(0, 1))
>>> x = sim.run(n_runs=1, method='nf')

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

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.
If obs_species_only is True, only the species contained within observables are output by cupSODA. All other concentrations are set to ‘nan’.

References

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.
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

An exception is thrown if tspan is not defined in either __init__`or `run.
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().
compiler: Choice of compiler for ODE system: cython, weave (Python 2 only), theano or python. Leave unspecified or equal to None for auto-select (tries weave, then cython, then python). Cython, weave and theano all compile the equation system into C code. Python is the slowest but most compatible.
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

In early versions of the Simulator class, tspan, initials and param_values supplied to this method persisted to future run() calls. This is no longer the case.

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.]

A_total trajectory for second run

>>> print(simulation_result.observables[1]['A_total'])         
[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

In early versions of the Simulator class, tspan, initials and param_values supplied to this method persisted to future run() calls. This is no longer the case.

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=None, observables_and_expressions=None, squeeze=True, simulations_per_param_set=1, model=None, initials=None, param_values=None)[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.
model: pysb.Model
initials: numpy.ndarray
param_values: numpy.ndarray: model, initials, param_values are an alternative constructor mechanism used when loading SimulationResults from files (see SimulationResult.load()). Setting just the simulator argument instead of these arguments is recommended.

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.

classmethod load(filename, dataset_name=None, group_name=None)[source]¶

Load a SimulationResult from a file (HDF5 format)

For a description of the file format see save()

Parameters:

filename: str: Filename from which to load data
dataset_name: str or None: Dataset name. Can be left as None when the group specified only contains one dataset, which will then be selected. If None and more than one dataset is in the group, a ValueError is raised.
group_name: str or None: Group name. This is typically the name of the model. Can be left as None when the file only contains one group, which will then be selected. If None and more than group is in the file a ValueError is raised.

Returns:

SimulationResult: Set of trajectories and associated metadata loaded from the file

nsims¶: The number of simulations in this SimulationResult

observable(pattern)[source]¶

Calculate a pattern’s trajectories without adding to model

This method calculates an observable “on demand” using any supplied MonomerPattern or ComplexPattern against the simulation result, without re-running the simulation.

Note that the monomers within the supplied pattern are reconciled with the SimulationResult’s internal copy of the model by name. This method only works on simulations which calculate species trajectories (i.e. it will not work on network-free simulations).

Raises a ValueError if the pattern does not match at least one species.

Parameters:	pattern: pysb.MonomerPattern or pysb.ComplexPattern An observable pattern to match
Returns:	pandas.Series Series containing the simulation trajectories for the specified observable

Examples

>>> from pysb import ANY
>>> from pysb.examples import earm_1_0
>>> from pysb.simulator import ScipyOdeSimulator
>>> simres = ScipyOdeSimulator(earm_1_0.model, tspan=range(5)).run()
>>> m = earm_1_0.model.monomers

Observable of bound Bid:

>>> simres.observable(m.Bid(b=ANY))
time
0    0.000000e+00
1    1.190933e-12
2    2.768582e-11
3    1.609716e-10
4    5.320530e-10
dtype: float64

Observable of AMito bound to mCytoC:

>>> simres.observable(m.AMito(b=1) % m.mCytoC(b=1))
time
0    0.000000e+00
1    1.477319e-77
2    1.669917e-71
3    5.076939e-69
4    1.157400e-66
dtype: float64

observables¶: List of trajectory sets. The first dimension contains observables.

save(filename, dataset_name=None, group_name=None, append=False, include_obs_exprs=False)[source]¶

Save a SimulationResult to a file (HDF5 format)

HDF5 is a hierarchical, binary storage format well suited to storing matrix-like data. Our implementation requires the h5py package.

Each SimulationResult is treated as an HDF5 dataset, stored within a group which is specific to a model. In this way, it is possible to save multiple SimulationResults for a specific model.

A group is first created in the HDF file root (see group_name argument). Within that group, a dataset “_model” has a pickled version of the PySB model. SimulationResult are stored as groups within the model group.

The file hierarchy under group_name/dataset_name/ then consists of the following HDF5 gzip compressed HDF5 datasets: trajectories, param_values, initials, tout, observables (optional) and expressions (optional); and the following attributes: simulator_class (pickled Class), simulator_kwargs (pickled dict), squeeze (bool), simulations_per_param_set (int), pysb_version (str), timestamp (ISO 8601 format).

Custom attributes can be stored in the SimulationResult’s custom_attrs dictionary. Keys should be strings, values can be any picklable object. When saved to HDF5, these custom attributes will be prefixed with usrattr_.

Parameters:

filename: str: Filename to which the data will be saved
dataset_name: str or None: Dataset name. If None, it will default to ‘result’. If the dataset_name already exists within the group, a ValueError is raised.
group_name: str or None: Group name. If None, will default to the name of the model.
append: bool: If False, raise IOError if the specified file already exists. If True, append to existing file (or create if it doesn’t exist).
include_obs_exprs: bool: Whether to save observables and expressions in the file or not. If they are not included, they can be recreated from the model and species trajectories when loaded back into PySB, but you may wish to include them for use with external software, or if you have complex expressions which take a long time to compute.

species¶: List of trajectory sets. The first dimension contains species.

class pysb.simulator.PopulationMap(complex_pattern, lumping_rate, counter_species=None)[source]¶

Population map for BioNetGen hybrid particle/population simulation

For use with the BngSimulator.

References

Hogg et al. 2014: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003544

BioNetGen HPP documentation: http://bionetgen.org/index.php/Hybrid_particle-population_model_generator

Simulation tools (pysb.simulator)¶

Simulation tools (`pysb.simulator`)¶