Introduction¶
PySCeS: the Python Simulator for Cellular Systems is an extendable toolkit for the analysis and investigation of cellular systems. It is available for download at http://pysces.sourceforge.net
Welcome! This users guide will get you started with the basics of modelling cellular systems with PySCeS. It is meant to be used together with the input file guide. If you already have PySCeS installed continue straight on if not, Installing and configuring contains instructions on building and installing PySCeS.
PySCeS is distributed under the PySCeS (BSD style) license and is made freely available as Open Source software. See LICENCE.txt for details.
PySCeS continued development depends, to a large degree on support and feedback from Systems Biology community, if you use PySCeS in your work please cite it using the following reference:
Brett G. Olivier, Johann M. Rohwer and Jan-Hendrik S. Hofmeyr Modelling cellular systems with PySCeS, Bioinformatics, 21, 560-561, DOI 10.1093/bioinformatics/bti046.
We hope that you will enjoy using our software. If, however, you find any unexpected features (i.e. bugs) or have any suggestions on how we can improve PySCeS please let us know.
The PySCeS development team.
Getting started¶
Loading PySCeS¶
In this section we assume you have PySCeS installed and
configured (see Installing and configuring for details) and a
correctly formatted PySCeS input file that describes a cellular
system in terms of its reactions, species and parameters. For a
detailed description of the PySCeS Model Description Language
see the PySCeS Input File Guide. Note that on all platforms
PySCeS model files have the extension .psc
.
To begin modelling we need to start up an interactive Python shell
(we suggest iPython) and load PySCeS with import pysces
:
Python 2.7.11 |Anaconda 2.4.1 (64-bit)| (default, Dec 7 2015, 14:10:42) [MSC v.1500 64 bit (AMD64)]
Type "copyright", "credits" or "license" for more information.
In [1]: import pysces
Matplotlib backend set to: "TkAgg"
Matplotlib interface loaded (pysces.plt.m)
Continuation routines available
NLEQ2 routines available
You are using NumPy (1.10.1) with SciPy (0.16.0)
RateChar is available
Parallel scanner is available
PySCeS environment
******************
pysces.model_dir = d:\\mypysces\\pscmodels
pysces.output_dir = c:\\Pysces
***********************************************************************
* Welcome to PySCeS (0.9.3) - Python Simulator for Cellular Systems *
* http://pysces.sourceforge.net *
* Copyright(C) B.G. Olivier, J.M. Rohwer, J.-H.S. Hofmeyr, 2004-2016 *
* Triple-J Group for Molecular Cell Physiology *
* Stellenbosch University, ZA and VU University Amsterdam, NL *
* PySCeS is distributed under the PySCeS (BSD style) licence, see *
* LICENCE.txt (supplied with this release) for details *
***********************************************************************
PySCeS is now ready to use. If you would like to test your installation try running the test suite:
pysces.test()
this also copies the test models supplied with PySCeS into your model directory.
Creating a PySCeS model object¶
This guide uses the test models supplied with PySCeS as examples, if you would like to use them and have not already done so run the PySCeS tests (described in the previous section).
Before modelling, a PySCeS model object needs to be instantiated.
As a convention we use mod
as the instantiated model
instance. The following code creates such an instance using the
test input file, pysces_test_linear1.psc
:
>>> mod = pysces.model('pysces_test_linear1')
Assuming extension is .psc
Using model directory: C:\mypysces\pscmodels
C:\mypysces\pscmodels\pysces_test_linear1.psc loading .....
Parsing file: C:\mypysces\pscmodels\pysces_test_linear1.psc
Calculating L matrix . . . . . . . done.
Calculating K matrix . . . . . . . done.
When instantiating a new model object, PySCeS input files are
assumed to have a .psc
extension. If the specified input
file does not exist in the input file directory (e.g.
misspelled filename) a list of existing input files is shown
and the user is given an opportunity to enter the correct
filename.
Advanced¶
The model constructor can also be used to specify a model directory other than the default model path:
mod = pysces.model('pysces_test_linear1', dir='c:\\Pysces\\psc')
alternatively input files can also be loaded from a string:
>>> F = file('c:\\Pysces\\psc\\pysces_test_linear1.psc', 'r')
>>> pscS = F.read()
>>> F.close()
>>> mod = pysces.model('test_lin1s', loader='string', fString=pscS)
Assuming extension is .psc
Using model directory: C:\mypysces\pscmodels
Using file: test_lin1s.psc
C:\mypysces\pscmodels\orca\test_lin1s.psc loading .....
note that now the input file is saved and loaded as
model_dir\orca\test_lin1s.psc
.
Creating a PySCeS model object¶
Once a new model object has been created it needs to be loaded. During the load process the input file is parsed, the model description is translated into Python data structures and a stoichiometric structural analysis is performed:
**New** in PySCeS 0.7.1+ model loading is now automatically performed
when the model object is instantiated. This behaviour is controlled
by the *autoload* argument (default = True). To keep backwards compatibility
with older modelling scripts the first time ``doLoad()`` is called a warning
is generated. Subsequent calls to ``doLoad()`` will reload the model as required.
Once loaded, all the model elements contained in the input file
are made available as model (mod
) attributes so that in the
input file where you might find initialisations such as s1 =
1.0
and k1 = 10.0
these are now available as mod.s1
and mod.k1
. For variable species and compartments an
additional attribute is created which contains the elements
initial (as opposed to current) value. These are constructed as
<name>*_init*:
>>> mod.s1
1.0
>>> mod.s1_init
1.0
>>> mod.k1
10.0
Any errors generated during the loading process (almost always) occur as a result of syntax errors in the input file. These error messages may not be intuitive for example, ‘list out of range’ exception usually indicates a missing multiplication operator “3(” instead of “3*(” or unbalanced parentheses.
Basic model attributes¶
Some basic model properties are accessible once the model is loaded:
mod.ModelFile
, the name of the model file that was used.mod.ModelDir
, the input file directory.mod.ModelOutput
, the PySCeS work/output directory.- Parameters are available as attributes directly as specified
in the input file e.g.
k1
ismod.k1
. - External (fixed) species are made available in the same way.
- Internal (variable) species are treated in a similar way except that an
additional attribute (parameter) is created to hold the species initial value
(as specified in the input file), e.g., from
s1
,mod.s1
andmod.s1_init
are instantiated as model object attributes. - Compartments are also are assigned an initial value.
- Rate equations are translated into objects that return their current value
when called
mod.R1()
.
All basic model attributes that are described here can be changed interactively. However, if the model rate equations need to be changed, this should be done in the input file after which the model should be re-instantiated and reloaded.
Groups of model properties (either tuples, lists or dictionaries)
mod.species
the model’s variable species names (ordered relative to the stoichiometric matrix rows).mod.reactions
reaction names and ordered to the stoichiometric matrices columns.mod.parameters
all parameters (including fixed species)mod.fixed_species
only the fixed species namesmod.__rate_rules__
a list of rate rules defined in the model
Advanced¶
These attributes are used by PySCeS to store additional information about the basic model components, generally they are supplied by the parser and should almost never be changed directly.
mod.__events__
a list of event object references which can be interrogated for event information. For example if you want a list of event names try[ev.name for ev in mod.__events__]
mod.__rules__
a dictionary containing information about all rules defined for this modelmod.__sDict__
a dictionary of species informationmod.__compartments__
a dictionary containing compartment information
Modelling¶
Structural Analysis¶
As part of the model loading procedure, doLoad()
automatically performs
a stoichiometric (structural) analysis of the model. The structural
properties of the model are captured in stoichiometric matrix (N),
kernel matrix (K) and link matrix (L). These matrices can
either be displayed with a mod.showX()
method or used in further
calculations as numeric arrays. The formal definition of these matrices,
as they are used in PySCeS, is described in [1].
The structural properties of a model are available in two forms, as new-style objects which have all the array properties neatly encapsulated or as legacy attributes. Although both exist it is highly recommended to use the new objects.
Structural Analysis - new objects¶
For alternate descriptions of these model properties see the next (legacy) section.
mod.Nmatrix
view withmod.showN()
mod.Nrmatrix
view withmod.showNr()
mod.Lmatrix
view withmod.showK()
mod.L0matrix
mod.Kmatrix
view withmod.showL()
mod.K0matrix
mod.showConserved()
displays any moiety conserved relationships (if present).mod.showFluxRelationships()
shows the relationships between dependent and independent fluxes at steady state.
All new structural objects have an array attribute which holds the actual NumPy array data as well as ridx and cidx which hold the row and column indexes (relative to the stoichiometric matrix) as well as the following methods:
.getLabels()
return the matrix labels as tuple([rows], [columns]).getColsByName()
extract column(s) with label.getRowsByName()
extract row(s) with label.getIndexes()
return the matrix indices (relative to the Stoichiometric matrix) as tuple((rows), (columns)).getColsByIdx()
extract column(s) referenced by index.getRowsByIdx()
extract row(s) referenced by index
Structural Analysis - legacy¶
mod.nmatrix
, N: displayed withmod.showN()
mod.kmatrix
, K: displayed withmod.showK()
mod.lmatrix
, L: displayed withmod.showL()
(an identity matrix means L does not exist i.e. no linear dependence).- If there are linear dependencies in the differential equations then the
reduced stoichiometric matrix of linearly independent, differential
equations Nr is available as
mod.nrmatrix
and is displayed withmod.showNr()
. If there is no dependence Nr = N. - In the case where there is linear dependence the moiety conservation sums
can be displayed by using
mod.showConserved()
. The conservation totals are calculated from the initial values of the variable species as defined in the model file. - When the K and L matrices exist, their dependent parts
(K0, L0) are available as
mod.kzeromatrix
andmod.lzeromatrix
. mod.showConserved()
prints any moiety conserved relationships (if present).mod.showFluxRelationships()
shows the relationships between dependent and independent fluxes at steady state.
If the mod.showX()
methods are used the row and column titles of the
various matrices are displayed with the matrix. Additionally, all of the
mod.showX()
methods accept an open file object as an argument. If this
file argument is present the method’s results are output to a file and not
printed to the screen. Alternatively, the order of each matrix dimension,
relative to the stoichiometric matrix, is available as either a row or
column array (e.g., mod.krow
, mod.lrow``, mod.kzerocol
).
Time simulation¶
PySCeS has interfaces to two ODE solvers either LSODA from ODEPACK (part of SciPy) or SUNDIALS CVODE (using PySundials). If PySundials is installed it will automatically select CVODE if compartments, events or rate rules are detected during model load as LSODA is not able capable of event handling or changing compartment sizes. If, however, you would like to select the solver manually this is also possible:
mod.mode_integrator = 'LSODA'
mod.mode_integrator = 'CVODE'
There are three ways of running a simulation:
Defining the start, end time and number of points and using the
mod.Simulate()
method directly:mod.sim_start = 0.0 mod.sim_end = 20 mod.sim_points = 50 mod.Simulate()
Using the
mod.doSim()
method where only the end time and points need to be specified. For example running a 20 point simulation from time 0 to 10:>>> mod.doSim(end=10.0, points=20.0)
Or using
mod.doSimPlot()
which runs the simulation and displays the results. In addition to doSim’s arguments the following arguments may be used:mod.doSimPlot(end=10.0, points=21, plot='species', fmt='lines', filename=None)
where:
- plot can be one of species, rates or all.
- fmt plot format, UPI backend dependent (default=’‘) or the CommonStyle ‘lines’ or ‘points’.
- filename if not None (default) then the plot is exported as filename.png
Another way of quickly visualising the results of a simulation
is to use the mod.SimPlot
method:
mod.SimPlot(plot='species', filename=None, title=None, log=None, format='lines')
where:
- plot: output to plot (default=’species’)
- ‘all’ rates and species
- ‘species’ species
- ‘rates’ reaction rates
- [‘S1’, ‘R1’, ] a list of model attributes (species, rates)
- filename (optional) if not None file is exported to filename (default=None)
- title the plot title (default=None)
- log use log axis for ‘x’, ‘y’, ‘xy’ (default=None)
- fmt plot format, UPI backend dependent (default=’‘) or the CommonStyle ‘lines’ or ‘points’.
Called without arguments, mod.SimPlot()
plots all the species
concentrations against time.
Simulation results¶
In PySCeS 0.7.x the simulation results have been consolidated
into a new mod.data_sim
object. By default species
concentrations/amounts, reaction rates and rate rules are
automatically added to the data_sim object. If extra
information (parameters, compartments, assignment rules) is
required this can easily be added using mod.CVODE_extra_output
, a
list containing any model attribute which is not added by default.
The mod.data_sim
object which has many methods for extracting simulation
data including:
data_sim.getTime()
return a vector of time pointsdata_sim.getSpecies()
returns array([[time], [species]])data_sim.getRates()
returns array([[time], [rates]])data_sim.getRules()
returns array([[time], [rate rules]])data_sim.getXData
returns array([[time], [CVODE_extra_output]])data_sim.getSimData(*args)
return an array consisting of time plus any available data series:mod.data_sim.getSimdata('s1', 'R1', 'Rule1', 'xData2')
data_sim.getAllSimData(*args)
return an array of all simulation datadata_sim.getDataAtTime(time)
return the results of the simulation at time.data_sim.getDataInTimeInterval(time, bound)
return the simulation data in the interval [time-bound, time+bound], if bound is not specified it is assumed to be the step size.
All the data_sim.get* methods by default only return a NumPy array containing the requested data, however if the argument lbls is set to True then both the array as well as a list of column labels is returned:
Sdata, Slabels = mod.data_sim.getSpecies(lbls=True)
This is very useful when using the PySCeS plotting interface (described later in this guide) to plot simulation results.
Advanced¶
PySCeS sets integrator options that attempt to configure the integration
algorithms to suit a particular model. However, almost every integrator
option can be overridden by the user.
Simulator settings are stored in PySCeS mod.__settings__
dictionary. For LSODA some useful keys are
(mod.__settings__[key]):
'lsoda_atol': 1.0e-012
'lsoda_rtol': 1.0e-007
'lsoda_mxordn': 12
'lsoda_mxords': 5
'lsoda_mxstep': 0
atol and rtol are the absolute and relative tolerances, while mxstep=0 means that LSODA chooses the number of steps (up to 500). If this is still not enough, PySCeS automatically increases the number of steps necessary to find a solution.
Additionally, CVODE allows per step step-size optimisation and automatic tolerance scaling:
'cvode_abstol': 1.0e-15
'cvode_abstol_factor': 1.0e-8
'cvode_auto_tol_adjust': True
'cvode_mxstep': 1000
'cvode_reltol': 1.0e-9
'cvode_stats': False
cvode_abstol is considered to be the minimum absolute tolerance, PySCeS first uses the initial species values multiplied by cvode_abstol_factor (so that [s]*[factor] >= [abstol]) to calculate its absolute tolerance. Once the simulation is underway PySCeS periodically readjusts the absolute tolerance on a per species basis based on the current species value.
If CVODE cannot find a solution in the given number of steps it automatically increases cvode_mxstep and tries again, however, it also keeps track of the number of times that this adjustment is required and if a specific threshold is passed it will begin to increase cvode_reltol by 1.0e3 (to a maximal value of 1.0e-3). Finally, if cvode_stats is enabled CVODE will display a report of its internal parameters after the simulation is complete.
Steady-state analysis¶
PySCeS solves for a steady state using either the non-linear solvers HYBRD, NLEQ2 or forward integration. By default PySCeS has solver fallback enabled which means that if a solver fails or returns an invalid result (i.e., contains negative concentrations) it switches to the next available solver. The solver chain is as follows:
- HYBRD (can handle ‘rough’ initial conditions, converges quickly).
- NLEQ2 (highly optimised for extremely non-linear systems, more sensitive to bad conditioning and slightly slower convergence).
- FINTSLV (finds a result when the change in max([species]) is less than 0.1%; slow convergence).
Solver fallback can be disabled by setting mod.mode_solver_fallback =
0
. Each of the three solvers is highly configurable and although the
default settings should work for most models configurable options
can be set in by way of the mod.__settings__ dictionary.
To calculate a steady state use the mod.doState()
method:
>>> mod.doState()
(hybrd) The solution converged.
The results of a steady-state evaluation are stored as arrays as well as
individual attributes and can be easily displayed using the
mod.showState()
method:
mod.showState()
displays the current steady-state values of both the species and fluxes.- For each reaction (e.g.
R2
) a new attributemod.J_R2
, which represents its steady-state value, is created. - Similarly, each species (e.g.
mod.s2
) has a steady-state attributemod.s2_ss
mod.state_species
inmod.species
order.mod.state_flux
inmod.reactions
order.
There are various ways of initialising the steady-state solvers although, in general, the default values can be used.
mod.mode_state_init
initialises the solver using either the initial values (0), a value close to zero (1). The default behaviour is to use the initial values.
New: mod.data_sstate¶
New to PySCeS 0.7 is the mod.data_sstate
object that by
default stores steady-state data (species, fluxes, rate rules)
in a manner similar to mod.data_sim. One notable exception is
that the current steady-state values are also made available as
attributes to this object (e.g. species S1’s steady-state value
is stored as mod.data_sstate.S1
). Using the
mod.STATE_extra_output
list it is possible to store user
defined data in the data_sstate object. Steady-state data can be
easily retrieved using the by now familiar .get* methods.
data_sstate.getSpecies()
returns a species arraydata_sstate.getFluxes()
returns a flux arraydata_sstate.getRules()
returns a rate rule arraydata_sstate.getXData()
returns an array defined in STATE_extra_outputdata_sstate.getStateData(*args)
return user defined array of data (‘S1’,’R2’)data_sstate.getAllStateData()
return all state data as an array
All these methods also accept the lbls=True argument in which case they return both array data and a label list:
ssdat, sslbl = mod.data_sstate.getSpecies(lbls=True)
Stability¶
PySCeS can analyse the stability of systems that can attain a steady state. It does this by calculating the Eigen values of the Jacobian matrix for the reduced system of independent ODE’s:
- ``mod.doEigen()`` calculates a steady-state and performs the stability analysis
- ``mod.showEigen`` prints out a stability report
- ``mod.doEigenShow()`` combines both of the above
The Eigen values are also available as attributes
mod.lambda1
etc. By default the Eigen values are stored as
mod.eigen_values
but if
mod.__settings__['mode_eigen_output'] = 1
is set both the
Eigen values as well as the left and right Eigen vectors are
stored as mod.eigen_vecleft
and mod.eigen_vecright
respectively. Please note that there is currently no guarantee
that the order of the Eigen value array corresponds to the
species order.
Metabolic Control Analysis¶
For practical purposes the following methods are collected into a set of meta-routines that all first solve for a steady state and then the required Metabolic Control Analysis (MCA) [2], [3] evaluation methods.
Elasticities¶
The elasticities towards both the variable species and parameters can be
calculated using mod.doElas()
which generates as output:
- Scaled elasticities generated as
mod.ecRate_Species
, e.g.mod.ecR4_s2
mod.showEvar()
displays the non-zero elasticities calculated with respect to the variable species.mod.showEpar()
displays the non-zero parameter elasticities.
As a prototype we also store the elasticities in an object,
mod.ec.*
this may become the default way of accessing
elasticity data in future releases but has not been stabilised
yet.
Control coefficients¶
Both control coefficients and elasticities can be calculated using a single
method, mod.doMca()
.
mod.showCC()
displays the complete set of flux and concentration control coefficients.- Individual control coefficients are generated as either
mod.ccSpecies_Rate
for a concentration control coefficient, e.g.mod.ccs1_R4
. - Similarly,
mod.ccJFlux_Rate
is a flux control coefficient e.g.mod.ccJR1_R4
.
As it is generally common practice to use scaled elasticities
and control coefficients PySCeS calculated these by default.
However, it is possible to generate unscaled elasticities and
control coefficients by setting the attribute
mod.__settings__['mode_mca_scaled'] = 0
in which case the
model attributes are attached as mod.uec
and mod.ucc
respectively.
As a prototype we also store the elasticities in an object,
mod.cc.*
this may become the default way of accessing
control coefficient data in future releases but has not been
stabilised yet.
Response coefficients¶
A new PySCeS feature is the ability to calculate the parameter response
coefficients for a model with the mod.doMcaRC()
method. Unlike the
elasticities and control coefficients the response coefficients are made
available as a single attribute mod.rc
. This attribute is a data
object, containing the response coefficients as attributes and has the
following methods:
rc.var_par
individual response coefficients can be accessed as attributes made up ofvariable_parameter
e.g.mod.rc.R1_k1
rc.get('var', 'par')
return a response coefficientrc.list()
returns all response coefficients as a dictionary of {key:value} pairsrc.select('attr', search='a')
select all response coefficients that refer to'attr'
e.g.select('R1')
orselect('k2')
rc.matrix
: the matrix of response coefficientsrc.row
: row labelsrc.col
: column labels
Responce coefficients with respect to moiety-conserved sums¶
The mod.doMcaRC()
method only calculates response coefficients with respect to explicit model parameters. However, in models with moiety-conservation the total concentration of all the species that form part of a particular moiety-conserved cycle is also a parameter of the model. PySCeS infers such moiety-conserved sums from the initial species concentrations specified by the user. In some cases it might be interesting to consider the effects that a change in the total concentration of a moiety will have on the steady-state. This analysis may be done with the method mod.doMcaRCT()
.
Since moiety-conserved sums are not explicitly named in PySCeS model files, 'T_'
is prepended to all the species names listed in mod.Consmatrix.row
. For instance, if the dependent species in a moiety-conserved cycle is 'A'
, then 'T_A'
designates the moiety-conserved sum.
The object mod.rc
is augmented with the results of mod.doMcaRCT()
. Response coefficients may thus be accessed with mod.rc.get('var', 'T_par')
.
Parameter scanning¶
Single dimension parameter scans¶
PySCeS has the ability to quickly generate and plot single dimension
parameter scans. Scanning a parameter typically involves changing a
parameter through a range of values and recalculating the steady state at
each step. Two methods are provided which simplify this task,
mod.Scan1()
is provided to generate the scan data while
mod.Scan1Plot()
is used to visualise the results. The first step is to
define the scan parameters:
mod.scan_in
is a string defining the parameter to be scanned e.g.'x0'
mod.scan_out
is a list of strings representing the attribute names you would like to track in the output eg.['J_R1','J_R2','s1_ss','s2_ss']
You also need to define the range of points that you would like to scan over. For a linear range SciPy has a useful function
scipy.linspace(start, end, points)
(SciPy can be accessed by typingimport scipy
in your Python shell). If you need to generate a log range usescipy.logspace(start, end, points)
.Both
scipy.linspace
andscipy.logspace
use the number of points (including the start and end points) in the interval as an input. Additionally, the start and end values ofscipy.logspace
must be entered as indices, e.g. to start the range at 0.1 and end it at 100 you would writescipy.logspace(-1, 2, steps)
. Setting up a PySCeS scan session might look something like:>>> import scipy >>> mod.scan_in = 'x0' >>> mod.scan_out = ['J_R1','J_R6','s2_ss','s7_ss'] >>> scan_range = scipy.linspace(0,100,11)
Before starting the parameter scan, it is important to check that all the
model attributes involved in the scan do actually exist. For example,
mod.J_R1
is created when mod.doState()
is executed, likewise all
the elasticities (mod.ecR_S
) and control coefficients (mod.ccJ_R
)
are only created when the mod.doMca()
method is called. If all the
attributes exist you can perform a parameter scan using the
mod.Scan1(scan_range)
method which takes your predefined scan range as
an argument:
>>> mod.Scan1(scan_range)
Scanning ...
11 (hybrd) The solution converged.
(hybrd) The solution converged ...
done.
When the scan has been successfully completed, the results are stored in
the array (mod.scan_res
) that has mod.scan_in
as its first column
followed by columns that represent the data defined in mod.scan_out
(if
invalid steady states are generated during the scan they are replaced by
NaN). Scan1 also reports the scan parameter values which generated the
invalid states.} . If one or more of your input or output parameters is not
a valid model attribute, it will be ignored. Once the parameter scan data
has been generated, the next step is to visualise it using the
mod.Scan1Plot()
method:
>>> mod.Scan1Plot(plot=[], title=None, log=None, format='lines', filename=None)
- plot if empty mod.scan_out is used, otherwise any subset of mod.scan_out (default=[])
- filename the filename of the PNG file (default=None, no export)
- title the plot title (default=None)
- log if None a linear axis is assumed otherwise one of [‘x’,’xy’,’xyz’] (default=None)
- format the backend dependent line format (default=’lines’) or the CommonStyle ‘lines’ or ‘points’.
Called without any arguments Scan1Plot plots all of mod.sim_out against mod.sim_in.
Two dimension parameter scans¶
Two dimension parameter scans can also easily be generated using the mod.Scan2D
method:
>>> mod.Scan2D(p1, p2, output, log=False)
- p1 is a list of [model parameter 1, start value, end value, points]
- p2 is a list of [model parameter 2, start value, end value, points]
- output the steady-state variable e.g. ‘J_R1’ or ‘A_ss’
- log if True scan using log ranges for both axes
To plot the results of two dimensional scan use the mod.Scan2DPlot
method. Note
that as Matplotlib cannot produce 3D plots the GnuPlot interface must be active
(see the section on plotting later on in this guide):
>>> mod.Scan2DPlot(title=None, log=None, format='lines', filename=None)
- filename the filename of the PNG file (default=None, no export)
- title the plot title (default=None)
- log if None a linear axis is assumed otherwise one of [‘x’,’xy’,’xyz’] (default=None)
- format the backend dependent line format (default=’lines’) or the CommonStyle ‘lines’ or ‘points’.
Multi-dimension parameter scans¶
This new PySCeS feature allows multi-dimensional parameter scanning. Any
combination of parameters is possible and can be added as master
parameters that change independently or slave parameters whose change is
coordinated with the previously defined parameter. Unlike mod.Scan1()
this function is accessed via the pysces.Scanner
class that is
instantiated with a loaded PySCeS model object:
>>> sc1 = pysces.Scanner(mod)
>>> sc1.addScanParameter('x3', 1, 10, 11)
>>> sc1.addScanParameter('k2', 0.1, 1000, 5, log=True)
>>> sc1.addScanParameter('k4', 0.1, 1000, 5, log=True, slave=True)
>>> sc1.addUserOutput('J_R1', 's1_ss')
>>> sc1.Run()
... scan: 55 states analysed
>>> sc1_res = sc1.getResultMatrix()
>>> print sc1_res[0]
array([1., 0.1, 0.1, 97.94286647, 49.1380999])
>>> print sc1_res[-1]
array([1.0e+01, 1.0e+03, 1.0e+03, -3.32564878e+00, 3.84227702e-03])
In this scan we define two independent (x3, k2
) and one dependent
(k3
) scan parameters and track the changes in the steady state
variables J_R1
and s1_ss
. Note that k2
and k4
use a
logarithmic scale. Once run the input parameters cannot be altered,
however, the output can be changed and the scan rerun.
sc1.addScanParameter(name, start, end, points, log, slave)
wherename
is the input parameter (as a string),start
andend
define the range with the required number ofpoints
. Whilelog
andslave
are boolean arguments indicating the point distribution and whether the axis is independent or not.sc1.addUserOutput(*args)
an arbitrary number of model attributes to be output can be added (this method automatically tries to determine the level of analysis necessary) e.g.addUserOutput('J_R1', 'ecR1_k2')
sc1.Run()
run the scan, if subsequent runs are required after changing output parameters usesc1.RunAgain()
. Note that it is not possible to change the input parameters once a scan has been run, if this is required a new Scanner object should be created.sc1.getResultMatrix(stst=False)
return the scan results as an array containing both input and output if stst = True append the steady-state fluxes and concentrations to the user output so that output has dimensions [scan_parameters]+[state_species+state_flux]+[Useroutput] otherwise return the default [scan_parameters]+[Useroutput].sc1.UserOutputList
the list of output namessc1.UserOutputResults
an array containing only the outputsc1.ScanSpace
the generated list of input parameters.
Plotting¶
The PySCeS plotting interface has been completely rewritten to facilitate the use of multiple plotting back-ends via a Unified Plotting Interface (UPI). Using the UPI we ensure that a specified subset of plotting methods is back-end independent (although the UPI can be extended with back-end specific methods). So far Matplotlib (default) and GnuPlot back-ends have been implemented.
The common UPI functionality is accessible as pysces.plt.*
while back-end specific functionality is available as
pysces.plt.m
(Matplotlib) and pysces.plt.g
(GnuPlot).
While the Matplotlib is activated by default GnuPlot needs to
be enabled (see Configuring PySCeS section) and then activated
using pysces.plt.p_activateInterface('gnuplot')
. All
installed interfaces can be activated or deactivated as
required:
>>> pysces.plt.p_activateInterface(interface)
>>> pysces.plt.p_deactivateInterface(interface)
where interface is either ‘matplotlib’ or ‘gnuplot’. The PySCeS UPI defines currently has the following methods:
plot(data, x, y, title='', format='')
plot a single line data[y] vs data[x]
- data the data array
- x x column index
- y y column index
- title is the line key
- format is the backend format string (default=’‘)
plotLines(data, x, y=[], titles=[], formats=[''])
plot multiple lines i.e. data[y1, y2, ] vs data[x]
- data the data array
- x x column index
- y is a list of line indexes, if empty all of y not including x is plotted
- titles a list of line keys, if empty Line1, Line2, etc is used
- formats a list (per line) of format strings, if formats only contains a single item, this format is used for all lines.
splot(data, x, y, z, title='', format='')
plot a surface i.e. data[z] vs data[y] vs data[x]
- data the data array
- x x column index
- y y column index
- z z column index
- title the surface key
- format a format string (default=’‘)
splotSurfaces(data, x, y, z=[], titles=[], formats=[''])
plot multiple surfaces i.e. data[z1, z2, ] vs data[y] vs data[x]
- data the data array
- x x column index
- y y column index
- z a list of z column indexes, if empty all data not including x, y are plotted
- titles a list of surface keys, if empty Surf1, Surf2 etc. is used
- formats is a list (per line) of format strings (default=’‘)
If formats only contains a single item, this format is used for all surfaces.
replot()
replot the current figure using all active interfaces (useful with GnuPlot type interfaces)
save(name, directory=None, dfmt='\%.8e')
save the plot data and (if possible) the back-end specific format file
- filename the filename
- directory optional (default = current working directory)
- dfmt the data format string (default=’%.8e’)
export(name, directory=None, type='png')
export the current plot as a <format> file (currently only PNG is guaranteed to be available on all back-ends).
- filename the filename
- directory optional (default = current working directory)
- type the file format (default=’png’).
setGraphTitle(title='PySCeS Plot')
set the graph title, unset if title=None
- title (string, default=’PySCeS Plot’) the graph title
setAxisLabel(axis, label='')
sets one or more axis label
- axis x, y, z, xy, xz, yz, zyx
- label label string (default=None)
Called with only the axis argument clears that axis’ label.
setKey(value=False)
enable or disable the current plot key, no arguments removes key.
- value boolean (default = False)
setLogScale(axis)
set axis to log scale
- axis is one of x, y, z, xy, xz, yz, zyx
setNoLogScale(axis)
set axis to a linear scale
- axis is one of x, y, z, xy, xz, yz, zyx
setRange(axis, min=None, max=None)
set one or more axis range
- axis is one of x, y, z, xy, xz, yz, zyx
- min is the range(s) lower bound (default=None, back-end auto-scales)
- max range(s) upper bound (default=None, back-end auto-scales)
setGrid(value)
enable or disable the graph grid
- value (boolean) True (on) or False (off)
plt.closeAll()
Close all active Matplolib figures
Displaying data¶
Displaying/saving model attributes¶
All of the showX()
methods, with the exception of mod.showModel()
operate in exactly the same way. If called without an argument, they
display the relevant information to the screen. Alternatively, if given an
open, writable (ASCII mode) file object as an argument, they write the
requested information to the open file. This allows the generation of
customised reports containing only information relevant to the model.
mod.showSpecies()
prints the current value of the model species (mod.M).mod.showSpeciesI()
prints the initial, parsed in, value of the model species (mod.Mi).mod.showPar()
prints the current value of the model parameters.mod.showState()
prints the current steady-state fluxes and species.mod.showConserved()
prints any moiety conserved relationships (if present).mod.showFluxRelationships()
shows the relationships between dependent and independent fluxes at steady statemod.showRateEq()
prints the reaction stoichiometry and rate equations.mod.showODE()
prints the differential equations.
Please note that the mod.showModel()
method is not
recommended for saving models as a PySCeS input file instead
use the Core2 based pysces.interface.writeMod2PSC
method
instead:
>>> pysces.interface.writeMod2PSC(mod, filename, directory, iValues=True, getstrbuf=False)
- filename: writes <filename>.psc or <model_name>.psc if None
- directory: (optional) an output directory
- iValues: if True (default) then the models initial values are used (or the current values if False).
- getstrbuf: if True a StringIO buffer is returned instead of writing to disk
Assuming you have loaded a model and run mod.doState()
the following
code opens a Python file object (rFile
), writes the steady-state
results to the file associated with the file object (results.txt
) and
then closes it again:
>>> rFile = file('results.txt','w')
>>> mod.showState() # print the results to screen
>>> mod.showState(rFile) # write the results to the file results.txt
>>> rFile.close()
Writing formatted arrays¶
The showX()
methods described in the previous sections allow the user a
convenient way to write the predefined matrices either to screen or file.
However, for maximum flexibility, PySCeS includes a suite of array writers
that enable one to easily write, in a variety of formats any array to a
file. Unlike the showX()
methods, the Write_array
methods are
specifically designed to write to data to a file.
In most modelling situations it is rare that an array needs to be stored or
displayed that does not have specific labels for its rows or columns.
Therefore, all the Write_array
methods take list arguments that can
contain either the row or column labels. Obviously, these lists should be
equal in length to the matrix dimension they describe and in the correct
order.
There are currently three custom array writing methods that work either with a 1D (vector) or 2D arrays (matrices). To allow an easy comparison of the output of these methods, all the following sections use the same example array as input.
Write_array()
¶
The basic array writer is the Write_array()
method. Using the default
settings this method writes a ‘tab delimited’ array to a file. It is
trivial to change this to a ‘comma delimited’ format by using the
separator = ' '
argument. Numbers in the array are formatted using the
global number format.
If column headings are supplied using the Col = []
argument they are
written above the relevant column and if necessary truncated to fit the
column width. If a column name is truncated it is marked with a *
and
the full length name is written as a comment after the array data.
Similarly row data can be supplied using the Row = []
argument in which
case the row names are displayed as a comment which is written after the
array data.
Finally, if the close_file
argument is enabled the supplied file object
is automatically closed after writing the array. The full call to the
method is:
>>> mod.Write_array(input, File=None, Row=None, Col=None, separator=' ')
which generates the array
## Write_array_linear1_11:12:23
#s0 s1 s2
-3.0043e-001 0.0000e+000 0.0000e+000
1.5022e+000 -5.0217e-001 0.0000e+000
0.0000e+000 1.5065e+000 -5.0650e-001
0.0000e+000 0.0000e+000 1.0130e+000
# Row: R1 R2 R3 R4
By default, each time an array is written, PySCeS includes an array header
consisting of the model name and the time the array was written. This
behaviour can be disabled by setting: mod.write_array_header = 0
Write_array_latex()
¶
The Write_array_latex()
method functions similarly to the generic
Write_array()
method except that it generates a formatted array that
can be included directly in a LaTeXdocument. Additionally, there is no
separator argument, column headings are not truncated and row labels appear
to the left of the matrix.
>>> mod.Write_array_latex(input, File=None, Row=None, Col=None)
which generates
%% Write_array_latex_linear1_11:45:03
\[
\begin{array}{r|rrr}
& $\small{s0}$ & $\small{s1}$ & $\small{s2}$ \\ \hline
$\small{R1}$ &-0.3004 & 0.0000 & 0.0000 \\
$\small{R2}$ & 1.5022 &-0.5022 & 0.0000 \\
$\small{R3}$ & 0.0000 & 1.5065 &-0.5065 \\
$\small{R4}$ & 0.0000 & 0.0000 & 1.0130 \\
\end{array}
\]
and in a typeset document appears as:
s0 s1 s2 R1 -0.3004 0.0000 0.0000 R2 1.5022 -0.5022 0.0000 R3 0.0000 1.5065 -0.5065 R4 0.0000 0.0000 1.0130
Installing and configuring¶
Before installing or building PySCeS the following software is required:
- Python 2.5 (or 2.4 plus the Elementree/cElementree packages)
- Numpy 1.2+
- SciPy 0.7.0 (0.6.x will work with NumPy > 1.0.5)
- GCC 4.2+ on Linux or MinGW GCC 3.4.5 on Windows is required for building PySCeS from source only
- Matplotlib 0.98.3 with the TkAgg backend (this is the default, but optional, plotting package but can be replaced with GnuPlot)
optional, but highly recommended, packages:
- libSBML 3.x install with the Python bindings for SBML support
- GnuPlot alternative plotting back-end
- iPython highly recommended for interactive modelling sessions
- SciTE editor for editing and running PySCeS based modelling programs
This software stack provides a powerful scientific programming platform which is used by PySCeS to provide a flexible Systems Biology Modelling environment.
PySCeS 0.7.0 itself has been modularised into a main package and a (growing) number of support modules which extends its core functionality. The most important of these is the advanced simulation support added by installing PySundials (http://pysundials.sf.net). Linux users should build and install the SUNDIALS library and PySundials (build instructions on the PySundials web site). Windows users can simply download and install the pysces_pysundials module.
- pysces_pysundials a binary port of SUNDIALS+PySundials for Windows
- pysces_metatool adds elementary mode support to PySCeS using MetaTool
- pysces_mariner SOAP based web services gateway, including a PySCeS server and remote client
- pysces_kraken (coming soon) PySCeS distributed processing module (currently distributed with PySCeS)
PySCeS and its extension modules use either the Python distutils or the Numpy distutils extensions. Assuming you have working versions of NumPy and SciPy on a Linux type operating systems building PySCeS is as easy as:
python setup.py install
On Windows (with MinGW) depending on your system configuration this becomes:
python setup.py config --compiler=mingw32 build --compiler=mingw32 install
In this release we have started prototyping Python egg support (currently only for windows) which is implemented via the setupegg.py build scripts.
By default PySCeS installs with a version of ZIB’s NLEQ2 non-linear solver. This software is distributed under its own non-commercial licence. Please see the README.txt document provided with this PySCeS installation for details.
Configuring¶
PySCeS has two configuration (*.ini) files that allows one to specify global (per installation) and local (per user options). Currently the multiuser options are only fully realised on Linux based systems. Global options are stored in the pyscfg.ini file which is created in your PySCeS installation directory (this is a Windows version with the Linux defaults indicated with in *value*):
[Pysces]
install_dir = c:\python25\lib\site-packages\pysces
gnuplot_dir = c:\model\gnuplot\binaries
model_dir = os.path.join(os.getenv('HOMEDRIVE')+os.path.sep,'Pysces','psc')
output_dir = os.path.join(os.getenv('HOMEDRIVE')+os.path.sep,'Pysces')
*model_dir = os.path.join(os.path.expanduser('~'),'Pysces','psc')*
*output_dir = os.path.join(os.path.expanduser('~'),'Pysces')*
The [Pysces] section contains information on the installation directory, the directory where the GnuPlot executable(s) can be found and the default model file and output directories. As we shall see some of these defaults can be overruled by the local configuration options:
[ExternalModules]
nleq2 = True
[PyscesModules]
pitcon = True
These sections define whether 3rd party algorithms (e.g. NLEQ2) are available for use, while the last section allows the alternate plotting backends to be enabled or disabled:
[PyscesConfig]
gnuplot = True
matplotlib = True
The user configuration files (pys_usercfg.ini) are created when
PySCeS is imported/run for the first time. On Windows this is
in <HOMEDRIVE>\Pysces
while on Linux this is in
$HOME\Pysces
. Once created the user configuration files can
be edited and will be used for every subsequent PySCeS session:
[Pysces]
output_dir = C:\mypysces
model_dir = C:\mypysces\pscmodels
gnuplot = False
Here I have customised my default model and output directories and disabled GnuPlot (enabled above). If required gnuplot_dir can also be set to point to an alternate location on a per user basis. Once you have PySCeS configured to your personal requirements you are ready to begin modelling.
References¶
Footnotes
[1] | Hofmeyr, J.-H.S. (2001) Metabolic control analysis in a nutshell, in T.-M. Yi, M. Hucka, M. Morohashi, and H. Kitano, eds, Proceedings of the 2nd International Conference on Systems Biology, pp. 291-300. |
[2] | Kacser, H. and Burns, J. A. (1973), The control of flux, Symp. Soc. Exp. Biol. 32, 65-104. |
[3] | Heinrich and Rappoport (1974), A linear steady-state treatment of enzymatic chains: General properties, control and effector strength, Eur. J. Biochem. 42, 89-95. |