PySCeS: the Python Simulator for Cellular
Systems is an extendable toolkit for the analysis and
investigation of cellular systems. PySCeS is available for download at
http://pysces.github.io and on GitHub where the source code is
maintained: https://github.com/PySCeS/pysces
Welcome! This user guide will get you started with the basics of modelling
cellular systems with PySCeS. It is meant to be read in conjunction with the
chapter on The PySCeS Model Description Language, which specifies the syntax of the input
file read by the program. 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) licence and is made
freely available as Open Source software. See LICENCE.txt for details.
The continued development of PySCeS 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 by opening an issue on Github.
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 (refer to
The PySCeS Model Description Language). 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 importpysces:
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 .pscUsing model directory: /home/jr/Pysces/psc/home/jr/Pysces/psc/pysces_test_linear1.psc loading .....Parsing file: /home/jr/Pysces/psc/pysces_test_linear1.pscCalculating 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.
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.
Note
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, whenever doLoad() is called a warning
is generated.
To force re-loading of a model from the input file, use mod.reLoad().
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 element’s
initial (as opposed to current) value. These are constructed as
<name>_init
>>> mod.s11.0>>> mod.s1_init1.0>>> mod.k110.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, 'listoutofrange' exception usually indicates a missing multiplication
operator(3( instead of 3*() or unbalanced parentheses.
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 is mod.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 and
mod.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, e.g. 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 ordered to the stoichiometric matrices
columns.
mod.parameters all parameters (including fixed species)
mod.fixed_species only the fixed species names
mod.__rate_rules__ a list of rate rules defined in the model
The following 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.nameforevinmod.__events__]
mod.__rules__ a dictionary containing information about all rules defined for this model
mod.__sDict__ a dictionary of species information
mod.__compartments__ a dictionary containing compartment information
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 the 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 NumPy 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.
For alternate descriptions of these model properties see the
next (legacy) section.
mod.Nmatrix view with mod.showN()
mod.Nrmatrix view with mod.showNr()
mod.Lmatrix view with mod.showL()
mod.L0matrix
mod.Kmatrix view with mod.showK()
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 indices (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
mod.lmatrix, L: displayed with mod.showL() (an identity
matrix means that no conservation relationships exist, i.e. there is no
linear dependence between species).
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 with
mod.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 and
mod.lzeromatrix.
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).
PySCeS has interfaces to two ODE solvers, either LSODA from
ODEPACK (part of SciPy) or SUNDIALS CVODE (using Assimulo).
If Assimulo is installed, PySCeS 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:
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)
Or using mod.doSimPlot() which runs the simulation and
graphically displays the
results. In addition to doSim()’s arguments the following arguments may
be used:
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.
Starting with PySCeS versions 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 has many methods for extracting simulation
data including:
data_sim.getTime() returns a vector of time points
data_sim.getAllSimData() return an array of all simulation data
data_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:
This is very useful when using the PySCeS plotting interface
(see Plotting) to plot simulation results.
For quick reference, simulation results are also available as a Numpy record
array (mod.sim). This allows the user to directly reference a particular
model attribute, e.g. mod.sim.Time, mod.sim.R1, or mod.sim.s1. Each
of these calls returns a vector of values of the particular model attribute
over the entire simulation (length of mod.sim_time).
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 the PySCeS mod.__settings__
dictionary. For LSODA some useful keys (default values indicated) are
(mod.__settings__[*key*]):
where 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.
The following are the most common options that can be set for CVODE,
with their defaults indicated:
where atol, rtol and mxstep are as above.
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). If cvode_stats is enabled CVODE will
display a report of its internal parameters after the
simulation is complete. Finally, CVODE will by default also output the time
points when events are triggered, even if these were not originally specified
in mod.sim_time. To disable this behaviour and strictly report only the
times in mod.sim_time, set cvode_return_event_timepoints to False.
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 (e.g., contains negative concentrations) it switches to the next
available solver. The solver chain is as follows:
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 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 attribute mod.J_R2, which
represents its steady-state value, is created.
Similarly, each species (e.g. mod.s2) has a steady-state attribute
mod.s2_ss.
mod.state_species is an array of steady-state species values in
mod.species order.
mod.state_flux is an array of steady-state fluxes in mod.reactions
order.
There are various ways of initialising the steady-state solvers although,
in general, the default values should be sufficient.
mod.mode_state_init initialises the solver using either the initial
values specified in the input file (0), or a value close to zero (1). The
default behaviour is to use the initial values.
Since PySCeS version 0.7 the mod.data_sstate object 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 array
data_sstate.getFluxes() returns a flux array
data_sstate.getRules() returns a rate rule array
data_sstate.getXData() returns an array defined in STATE_extra_output
data_sstate.getStateData(*args) return user defined array of data
('S1','R2')
data_sstate.getAllStateData() return all steady-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:
PySCeS can analyse the stability of systems that can attain a steady state.
It does this by calculating the eigenvalues of the Jacobian matrix for the
reduced system of independent ODEs.
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 eigenvalues are also available as attributes
mod.lambda1 etc. By default the eigenvalues are stored as
mod.eigen_values but if
mod.__settings__['mode_eigen_output']=1 is set, in addition to the
eigenvalues the left and right eigenvectors are
stored as mod.eigen_vecleft and mod.eigen_vecright,
respectively. Please note that there is currently no guarantee
that the order of the eigenvalue array corresponds to the
species order.
For ease of use the following methods are collected into a set of
meta-routines that all first solve for a steady state and then
perform the required
Metabolic Control Analysis (MCA) [2], [3] evaluation methods.
The elasticities towards both the variable species and parameters can be
calculated using mod.doElas() which generates as output:
Scaled elasticities referenced 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 fully stabilised
yet.
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 concentration-control coefficients are referenced as
mod.ccSpecies_Rate, 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 calculate 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 control coefficients in an object,
mod.cc.*; this may become the default way of accessing
control coefficient data in future releases but has not been fully
stabilised yet.
PySCeS can 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 of variable_parameter e.g. mod.rc.R1_k1
rc.get('var','par') return a response coefficient
rc.list() returns all response coefficients as a dictionary of
{key: value} pairs
rc.select('attr',search='a') select all response coefficients that
refer to 'attr' e.g. select('R1') or select('k2')
rc.matrix: the matrix of response coefficients
rc.row: row labels
rc.col: column labels
Response 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').
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.
'k0'
mod.scan_out is a list of strings representing the attribute names
to be tracked in the output, e.g.
['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 NumPy has a useful function
numpy.linspace(start,end,points) (NumPy can be accessed by importing it
in your Python shell via importnumpy). If you need to generate a log
range use numpy.logspace(start,end,points).
Both numpy.linspace and numpy.logspace use the number of points
(including the start and end points) in the interval as an input.
Additionally, the start and end values of numpy.logspace must be
entered as indices, e.g. to start the range at 0.1 and end it at 100 you
would write numpy.logspace(-1,2,steps). Setting up a PySCeS scan
session might look something like:
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 the specified input or output parameters are
not valid model attributes, they will be ignored. Once the parameter scan data
has been generated, the next step is to visualise it using the
mod.Scan1Plot() method:
plot if empty, mod.scan_out is used, otherwise any subset of mod.scan_out
(default= [])
filename the filename of the PNG file to save (default= None, no export)
title the plot title (default= None)
log if None a linear axis is assumed, otherwise one of
['x','y','xy'] (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.scan_out against
mod.scan_in.
In a similar way that simulation results are captured in the mod.sim array,
1D-scan results are also available as a Numpy record array (mod.scan) for
quick reference and easy access by the user. All the model attributes defined
in mod.scan_in and mod.scan_out can be accessed in this way, e.g.
mod.scan.x0, mod.scan.J_R1, mod.scan.s2_ss, etc.
Two-dimensional parameter scans can also easily be generated using the mod.Scan2D
method:
>>> mod.Scan2D(p1,p2,output,log=False)
p1 is a list of [modelparameter1,startvalue,endvalue,points]
p2 is a list of [modelparameter2,startvalue,endvalue,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: the GnuPlot interface must be active for this to work
(see the section on Plotting later on in this guide).
This PySCeS feature allows multi-dimensional parameter scanning. Any
combination of parameters is possible and can be added as leader
parameters that change independently or follower parameters whose change is
coordinated with the previously defined parameter. Unlike mod.Scan1()
this function is accessed via the pysces.Scanner class that is separately
instantiated with a loaded PySCeS model object:
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,follower) where
name is the input parameter (as a string), start and end define
the range with the required number of points, While log and
follower 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 attributes, use sc1.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 names
sc1.UserOutputResults an array containing only the output
sc1.ScanSpace the generated list of input parameters.
When performing large multi-dimensional parameter scans, PySCeS has the option
to perform the computation in parallel, either on a single machine with a
multi-core CPU, or on a multi-node cluster. This requires a working
ipyparallel installation (see also Installation). The functionality is accessed
via the pysces.ParScanner class, which has the same methods as the pysces.Scanner
class (see above) with a few multiprocessing-specific additions.
The parallel scanner class is instantiated with a loaded PySCeS model object:
>>> sc1=pysces.ParScanner(mod,engine='multiproc')
The additional engine argument specifies the parallel computation engine to
use:
'multiproc' - use Python’s internal multiprocessing module (default)
'ipcluster' - use ipcluster (refer to ipyparallel documentation)
There are two ways to run the scan:
sc1.Run() - runs the scan with a load-balancing task client; tasks are
queued and sent to nodes as these become available.
sc1.RunScatter() - compute tasks are evenly distributed amongst compute
nodes (“scattered”) and the results are returned (“gathered”) once all
the computations are complete. No load balancing is performed. May be
slightly faster than sc1.Run() if the individual tasks are very similar.
Not available withmultiproc!
Further input and output processing is as for pysces.Scanner. A few example
scripts illustrating the parallel scanning procedure are provided in the
pysces/examples folder of the installation.
The PySCeS plotting interface has written 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 Configuration section) and then activated
using pysces.plt.p_activateInterface('gnuplot'). All
installed interfaces can be activated or deactivated as
required:
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 values of the model species
(mod.M).
mod.showSpeciesI() prints the initial values of the model
species (mod.Mi), as parsed from the input file.
mod.showPar() prints the current values 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 state.
mod.showRateEq() prints the reaction stoichiometry and rate equations.
mod.showODE() prints the ordinary differential equations.
Note
The mod.showModel() method is not
recommended for saving models as a PySCeS input file,
use the Core2 based pysces.interface.writeMod2PSC method
instead:
filename: writes <filename>.psc or <model_name>.psc if None
directory: (optional) an output directory
iValues: if True (default) then the model initial values are used
(or the current values if False)
getstrbuf: if True a StringIO buffer is returned instead of writing to disk
For example, 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=open('results.txt','w')>>> mod.showState()# print the results to screen>>> mod.showState(rFile)# write the results to the file results.txt>>> rFile.close()
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 (matrix) array. To allow an easy comparison of
the output of these methods, all the following sections use the same
example array as input.
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:
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
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 LaTeX document. Additionally, there is no
separator argument, column headings are not truncated and row labels appear
to the left of the matrix.
PySCeS is developed primarily in Python and has been designed to
operate on multiple operating systems, i.e. Linux, Microsoft Windows and macOS.
PySCeS makes use of NumPy and SciPy for a number of functions
and needs a working SciPy stack (https://www.scipy.org)
to install and run.
IPython or the Jupyter notebook (optional, highly recommended for
interactive modelling sessions)
libSBML (optional). Python bindings for SBML support can be installed via
$ pip install python-libsbml
This software stack provides a powerful scientific programming
platform which is used by PySCeS to provide a flexible Systems
Biology Modelling environment.
PySCeS itself has been modularised into a main package
and a (growing) number of support modules which extend its
core functionality. It is highly recommended that
the following packages/modules are also installed:
Assimulo to enable CVODE support. This can be installed on Anaconda via the
conda-forge channel, or compiled from source
(https://jmodelica.org/assimulo).
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 https://github.com/PySCeS/pysces
for details.
Binary install packages for all three OSs and Python versions 3.7-3.10 are
provided. Anaconda users can conveniently install PySCeS with:
$ conda install -c conda-forge -c pysces pysces
Any dependencies will be installed automatically, including the optional dependencies
Assimulo, ipyparallel and libSBML.
Alternatively, you can use pip to install PySCeS from PyPI. Core
dependencies will be installed automatically.
$ pip install pysces
To install the optional dependences:
pipinstall"pysces[parscan]" - for ipyparallel
pipinstall"pysces[sbml]" - for libSBML
pipinstall"pysces[cvode]" - for Assimulo
pipinstall"pysces[all]" - for all of the above
Note
Installation of Assimulo via pip may well require C and Fortran compilers to be
properly set up on your system, as binary packages are only provided for a
very limited number of Python versions and operating systems on PyPI.
This is not guaranteed to work! In addition, the Assimulo version on PyPI is
severely outdated. If you require Assimulo, the conda
install is by far the easier option as up-to-date binaries are supplied
for all OS and recent Python versions.
All modern Linux distributions ship with gcc and gfortran. In addition, the
Python development headers (python-dev or python-devel, depending on your
distro) need to be installed.
Clone the source from Github as described above, change into the source
directory and run:
PySCeS has two configuration (*.ini) files that allow one to
specify global (per installation) and local (per user) options.
Global options are stored in the
pyscfg.ini file which is created in your PySCeS
installation directory upon install. The example below is a Windows
version; the exact values of install_dir and gnuplot_dir (if available)
will depend on your individual OS and Python setup and are determined on
install.
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.
This section also contains two further key-value pairs. If silentstart
(default False) is set to True, informational messages about the PySCeS
installation are not printed to the console on startup. The key
change_dir_on_start specifies if the working directory should be changed to
the PySCeS output directory (typically $HOME/Pysces or
%USERPROFILE%\Pysces) on startup. When set to False (the default), the
working directory is not changed.
As we shall see some of these defaults can be overridden by the local
configuration options.
These sections define whether third-party algorithms (NLEQ2 and PITCON)
are available for use, while the last section allows the alternate
plotting backends to be enabled or disabled:
[PyscesConfig]gnuplot=Truematplotlib=True
The user configuration file (pys_usercfg.ini) is created when
PySCeS is imported/run for the first time. On Windows this is
in %USERPROFILE%\Pysces while on Linux and macOS this is in
$HOME/Pysces. Once created, the user configuration files can
be edited and will be used for every subsequent PySCeS session.
For example, the above user configuration on a Windows system customises the
default model and output directories
and disables GnuPlot (enabled globally above). If required, gnuplot_dir
can also be set to point to an alternate location on a per-user
basis. The configuration keys silentstart and change_dir_on_start can also
be overridden here on a per-user basis.
Once you have PySCeS configured to your personal
requirements you are ready to begin modelling.
PySCeS: the Python Simulator for Cellular
Systems is an extendable toolkit for the analysis and investigation of
cellular systems. It is available for download from: http://pysces.github.io.
This section deals with the PySCeS Model Description Language (MDL), a
human-readable format for defining kinetic models.
PySCeS uses an ASCII text based input file to describe a
cellular system in terms of its stoichiometry, kinetics,
compartments and parameters. Input files may have any filename
with the single restriction that, for cross platform
compatibility, they must end with the extension .psc. Since version 0.7,
the PySCeS MDL has been updated and extended to be compatible with models
defined in the SBML standard.
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 by opening an issue on Github.
The basic description of a kinetic model in the PySCeS MDL contains
the following information:
whether any fixed (boundary) species are present
the reaction network stoichiometry
rate equations for each reaction step
parameter and boundary species initial values
the initial values of the variable species
Although it is in principle possible to define an ODE based model
without reactions or free species, for practical purposes PySCeS
requires a minimum of a single reaction. Once this information is
obtained it can be organised and written as a PySCeS input file.
While this list is the minimum information required for a PySCeS input
file, the MDL allows the definition of advanced models that contain
compartments, global units, functions, rate and assignment rules.
Since PySCeS 0.7 it is possible to define keywords that
specify model information. Keywords have the general form
<keyword>:<value>
The Modelname (optional) keyword, containing only
alphanumeric characters (or _), describes the model filename
(typically used when the model is exported via the PySCeS
interface module) while the Description keyword is a (short)
single line model description.
PySCeS 0.7+ supports the (optional) definition of a set of
global units. In doing so we have chosen to follow the general
approach used in the Systems Biology Modelling Language (SBML
L2V3) specification. The general definition of a PySCeS unit
is: `<UnitType>:<kind>,<multiplier>,<scale>,<exponent>`
where kind is a string describing the base unit (for SBML
compatibility this should be an SI unit) e.g. mole, litre,
second or metre. The base unit is modified by the multiplier,
scale and index using the following relationship:
<multiplier> * (<kind> * 10**<scale>)**<index>. The
default unit definitions are:
Symbol names (i.e. reaction, species, compartment, function,
rule and parameter names etc.) must start with either an
underscore or letter and be followed by any combination of
alpha-numeric characters or an underscore. Like all other
elements of the input file names are case sensitive:
R1_subApar1bext_1
Explicit access to the “current” time in a time simulation is
provided by the special symbol _TIME_. This is useful in
the definition of events and rules (see section on Advanced
model construction for more details).
Comments can be placed anywhere in the input file in one of two
ways, as single line comment starting with a # or as a Python-styled
multi-line triple quoted “””<comment>”””:
# everything after this is ignored"""This is a commentspread over afew lines."""
By default (as is the case in all PySCeS versions < 0.7) PySCeS
assumes that the model exists in a single unit volume
compartment. In this case it is not necessary to define a
compartment and the ODEs therefore describe changes in
concentration per time. However, if a compartment is defined,
PySCeS assumes that the ODEs describe changes in substance amount per
time in keeping with the SBML standard. Doing this affects
how the model is defined in the input
file (especially with respect to the definitions of rate
equations and species) and the user is strongly advised to
read the User Guide before building models in this way. The
general compartment definition is Compartment:<name>,<size>,<dimensions>, where <name> is the unique
compartment id, <size> is the size of the compartment (i.e.
length, volume or area) defined by the number of <dimensions>
(e.g. 1,2,3):
A relatively recent addition to the PySCeS MDL is the ability to define
SBML-styled functions. Simply put these are code substitutions that
can be used in rate equation definitions to, for example,
simplify the kinetic law. The general syntax for a function is
Function:<name>,<arglist>{<formula>} where <name> is the
unique function id, <arglist> is one or more comma separated
function arguments. The <formula> field, enclosed in curly
braces, may only make use of arguments listed in the
<arglist> and therefore cannot reference model attributes
directly. If this functionality is required a forcing function/assignment rule
(see Assignment rules) may be what you are looking
for.
Boundary species, also known as fixed or external species, are
a special class of parameter used when modelling biological
systems. The PySCeS MDL fixed species are declared on a single
line as FIX:<fixedlist>. The <fixedlist> is a space-separated list of
symbol names which should be initialised like
any other species or parameter:
FIX:Fru_exGlc_exATPADPUDPphosglycolysisSuc_vac
If no fixed species are present in the model then this
declaration should be omitted entirely.
The reaction stoichiometry and rate equation are defined together
as a single reaction step. Each step in the system is defined as
having a name (identifier), a stoichiometry (substrates are
converted to products) and rate equation (the catalytic activity,
described as a function of species and parameters). All reaction
definitions should be separated by an empty line. The general format
of a reaction in a model with no compartments is:
<name>:<stoichiometry><rateequation>
The <name> argument follows the syntax as discussed in
Symbol names and comments above;
however, when more than one compartment has
been defined it is important to locate the reaction in its
specific compartment. This is done using the @ operator:
where <compartment> is a valid compartment name. In either
case this then followed (either directly or on the next line)
by the reaction stoichiometry.
Each <stoichiometry> argument is defined in terms of reaction
substrates, appearing on the left hand side, and products on the
right hand side of an identifier which labels the reaction as
either reversible (=) or irreversible (>). If required each
reagent’s stoichiometric coefficient (PySCeS accepts both
integer and floating point) should be included in curly braces
{} immediately preceding the reagent name. If these are
omitted a coefficient of one is assumed.
{2.0}Hex_P=Suc6P+UDP# reversible reaction
Fru_ex>Fru# irreversible reaction
species_5>$pool# a reaction to a sink
The PySCeS MDL also allows the use of the $pool token that
represents a placeholder reagent for reactions that have no
net substrate or product. The reversibility of a reaction is only
used when exporting the model to other formats (such as SBML)
and in the calculation of elementary modes. It does not affect
the numerical evaluation of the rate equations in any way.
Central to any reaction definition is the <rate equation>
(SBML kinetic law). This should be written as valid Python
expression and may fall across more than one line. Standard
Python operators +-*/** are supported (note the Python
power e.g. 2^4 is written as 2**4). There is no shorthand
for multiplication with a bracket so -2(a+b)^h would be written as
-2*(a+b)**h and normal operator precedence applies:
+, -
addition, subtraction
*, /
multiplication, division
+x,-x
positive, negative
**
exponentiation
Operator precedence increases from top to bottom and left to
right (adapted from the Python Reference Manual).
The PySCeS MDL parser has been developed to parse and translate different
styles of infix into Python/Numpy-based expressions. The following
functions are supported in any mathematical expression:
log, log10, ln, abs (note, log is defined as natural
logarithm, equivalent to ln)
pow, exp, root, sqrt
sin, cos, tan, sinh, cosh, tanh
arccos, arccosh, arcsin, arcsinh, arctan, arctanh
floor, ceil, ceiling, piecewise
notanumber, pi, infinity, exponentiale
Logical operators are supported in rules, events, etc., but not
in rate equation definitions. The PySCeS parser understands
Python infix as well as libSBML and NumPy prefix notation.
andorxornot
>gt(x,y)greater(x,y)
<lt(x,y)less(x,y)
>=ge(x,y)geq(x,y)greater_equal(x,y)
<=le(x,y)leq(x,y)less_equal(x,y)
==eq(x,y)equal(x,y)
!=neq(x,y)not_equal(x,y)
Note that currently the MathML delay and factorial functions
are not supported. Delay is handled by simply removing it from
any expression, e.g. delay(f(x), delay) would be parsed as
f(x). Support for piecewise has been recently added
to PySCeS and will be discussed in the advanced features section (see
Piecewise).
A reaction definition when no compartments are defined:
When compartments are defined, note how now the reaction is now
given a location. This is because the ODEs formed from these
reactions must be in changes in substance (amount) per time, thus the rate
equation is multiplied by its compartment size. In this
particular example the species symbols represent concentrations
(Species_In_Conc: True):
Please note that at this time we are not certain if this form
of rate equation is translatable into valid SBML in a way that is
interoperable with other software.
The general form of any species (fixed, free) and parameter is
simply:
property=value
Initialisations can be written in any order anywhere in the
input file but for enhanced human readability these are usually
placed after the reaction that uses them or grouped at the end
of the input file. Both decimal and scientific notation are
allowed with the following provisions that neither floating
point (1.) nor scientific shorthand (1.e-3) syntax should
be used, instead use the full form (1.0e-3, 0.001 or
1.0).
Variable or free species are initialised differently depending
on whether compartments are present in the model or not. Although the
variable species concentrations are determined by
the parameters of the system, their initial values are used in
various places, calculating total moiety concentrations (if
present), time simulation initial values (e.g. time=zero) and
as initial guesses for the steady-state algorithms. If an empty
initial species pool is required it is not recommended to
initialise these values to zero (in order to prevent potential
divide-by-zero errors) but rather to a small value (e.g.
1.0e-8).
For a model with no compartments these initial values are assumed
to be concentrations:
NADH=0.001ATP=2.3e-3sucrose=1
In a model with compartments it is expected that the species
are located in a compartment (even if Species_In_Conc: False);
this is done using the @ symbol:
s1@Memb=0.01s2@Cell=2.0e-4
A word of warning, the user is responsible for making sure that
the units of the initialised species match those of the model.
Please keep in mind that all species (and anything that
depends on them) are defined in terms of the Species_In_Conc
keyword. For example, if the preceding initialisations were for
R1 (see Reaction section) then they would be concentrations
(as Species_In_Conc: True). However, in the next example, we
are initialising species for R4 and they are therefore in
amounts (Species_In_Conc: False).
s3@Memb=1.0s4@Cell=2.0
Fixed species are defined in a similar way and although they
are technically parameters, they should be given a location in
compartmental models:
# InitExtX0=10.0X4@Cell=1.0
However, fixed species are true parameters in the sense that
their associated compartment size does not affect their value
when it changes size. If compartment size-dependent behaviour
is required, an assignment or rate rule should be considered.
Finally, the parameters should be initialised. PySCeS checks if
a parameter is defined that is not present in the rate
equations and if such parameter initialisations are detected a
harmless warning is generated. If, on the other hand, an
uninitialised parameter is detected a warning is generated and
a value of 1.0 assigned:
Assignment rules or forcing functions are used to set the value
of a model attribute before the ODEs are evaluated. This model
attribute can either be a parameter used in the rate equations
(this is traditionally used to describe an equilibrium block), a
compartment, or an arbitrary parameter (commonly used to define
some sort of tracking function). Assignment rules can access
other model attributes directly and have the generic form !F<par>=<formula>, where <par> is the parameter assigned
the result of <formula>. Assignment rules can be defined
anywhere in the input file:
These rules would set the value of <par>, whose value
can be followed using the simulation and steady state
extra_output functionality (see Simulation results and
The steady-state data object).
PySCeS includes support for rate rules, which are
essentially directly encoded ODEs that are evaluated after
the ODEs defined by the model stoichiometry and rate
equations. Unlike the SBML rate rule, PySCeS allows one to
directly access a reaction symbol in the rate rules (this is
automatically expanded when the model is exported to SBML). The
general form of a rate rule is RateRule:<name>{<formula>}, where <name> is the model attribute (e.g.
compartment or parameter) whose rate of change is described by
the <formula>. It may also be defined anywhere in the input
file:
Time-dependent events may be defined. Since PySCeS 1.1.0 their definition
follows the event framework described in the SBML L3V2
specification. The general form of an event is:
As can be seen, an event consists of essentially two parts, a conditional
<trigger>, and a set of one or more <assignments>. Assignments have the
general form <par>=<formula>. Events have access to the “current”
simulation time using the _TIME_ symbol:
In order for PySCeS to handle events it is necessary to
have Assimulo installed (refer to General requirements).
In line with the SBML L3V2 specification, three optional keyword arguments
can be defined as a comma-separated list following the trigger definition. The
general syntax is <attribute>=<value>. The keywords are:
delay (float): specifies a delay between
when the trigger is fired (and the assignments are evaluated)
and the eventual assignment to the model. If this keyword is not specified, a
value of 0.0 is assumed.
priority (integer or None): specifies a priority for events that trigger at
the same simulation time. Events with a higher priority are executed before
those with a lower priority, while events without a priority (None) are
executed in random positions in the sequence. If this keyword is not
specified, a value of None is assumed.
persistent (boolean): is only relevant to events with a delay, where the
situation may occur that the trigger condition no
longer holds by the time the delay in the simulation has passed. The
persistent attribute specifies how to deal with this situation: if
True, the event executes nevertheless; if False, the event does not
execute if the trigger condition is no longer valid. If the keyword is not
specified, a default of True is assumed.
The following event illustrates the use of a delay of ten time units with a
non-persistent trigger and a priority of 3. In addition, the prefix notation
(used by libSBML) for the trigger is illustrated (PySCeS understands both
notations):
The legacy event specification (PySCeS versions <1.1), which did not
include keywords for priority and persistent, and in which the
delay was specified as a third positional argument (without
keyword), is still supported but deprecated:
Event:<name>,<trigger>,<delay>{<assignments>}
The following SBML event attributes are not implemented:
event.use_values_from_trigger_time=False For an event with a
delay, PySCeS always uses the assignment values from the time
when the event is triggered. When loading an SBML model
(see SBML import and export) that uses the
assignment values from the time of event firing and thus has this
event attribute set, a NotImplementedError is raised.
trigger.initial_value=False PySCeS always assumes that the initial
value of a trigger is True, i.e. the event cannot fire at time zero,
but that the simulation has to run for at least one iteration before the
trigger can be fired. When loading and SBML model that has the initial
value set to False, a NotImplementedError is raised.
Although technically an operator, piecewise functions are
sufficiently complicated to warrant their own section. A
piecewise operator is essentially an if, elif, …, else
logical operator that can be used to conditionally “set” the
value of some model attribute. Currently piecewise is supported
in rule constructs and has not been tested directly in rate
equation definitions. The piecewise function’s most basic
incarnation is piecewise(<val1>,<cond>,<val2>), which is
evaluated as:
In principle there is no limit on the number of conditional
statements present in a piecewise function; the condition can
be a compound statement (aorbandc) and may include the
_TIME_ symbol.
Some models contain reactions that are defined as only having substrates or
products, with the fixed (external) species not specified:
R1:A+B>R2:>C+D
The implication is that the relevant reagents appear from or disappear into
a constant pool. Unfortunately the PySCeS parser does not accept
such an unbalanced reaction definition and requires these pools to be
represented with a $pool token:
R1: A + B > $pool
R2: $pool > C + D
$pool is neither counted as a reagent nor does it ever appear in the
stoichiometry (think of it as dev/null) and no other $<str> tokens are
allowed.
SBML models can be imported into PySCeS by first converting them to the input file
(*.psc) format and then loading the input file into PySCeS as usual.
The conversion is done with the pysces.interface.convertSBML2PSC() method:
PySCeS test model: pysces_test_linear1.psc - this file is distributed with
PySCeS and copied to your model directory (typically $HOME/Pysces/psc) after
installation, when running pysces.test() for the first time.
This model includes the use of Compartments, KeyWords,
Units and Rules:
Modelname: MWC_wholecell2c
Description: Surovtsev whole cell model using J-HS Hofmeyr's framework
Species_In_Conc: True
Output_In_Conc: True
# Global unit definition
UnitVolume: litre, 1.0, -3, 1
UnitSubstance: mole, 1.0, -6, 1
UnitTime: second, 60, 0, 1
# Compartment definition
Compartment: Vcyt, 1.0, 3
Compartment: Vout, 1.0, 3
Compartment: Mem_Area, 5.15898, 2
FIX: N
R1@Mem_Area: N = M
Mem_Area*k1*(Pmem)*(N/Vout)
R2@Vcyt: {244}M = P # m1
Vcyt*k2*(M)
R3@Vcyt: {42}M = L # m2
Vcyt*k3*(M)*(P)**2
R4@Mem_Area: P = Pmem
Mem_Area*k4*(P)
R5@Mem_Area: L = Lmem
Mem_Area*k5*(L)
# Rate rule definition
RateRule: Vcyt {(1.0/Co)*(R1()+(1-m1)*R2()+(1-m2)*R3()-R4()-R5())}
RateRule: Mem_Area {(sigma_P)*R4() + (sigma_L)*R5()}
# Rate rule initialisation
Co = 3.07e5 # uM p_env/(R*T)
m1 = 244
m2 = 42
sigma_P = 0.00069714285714285711
sigma_L = 0.00012
# Assignment rule definition
!F S_V_Ratio = Mem_Area/Vcyt
!F Mconc = (M)/M_init
!F Lconc = (L)/L_init
!F Pconc = (P)/P_init
# Assignment rule initialisations
M_init = 199693.0
L_init = 102004
P_init = 5303
Mconc = 1.0
Lconc = 1.0
Pconc = 1.0
# Species initialisations
N@Vout = 3.07e5
Pmem@Mem_Area = 37.38415
Lmem@Mem_Area = 8291.2350678770199
M@Vcyt = 199693.0
L@Vcyt = 102004
P@Vcyt = 5303
# Parameter initialisations
k1 = 0.00089709
k2 = 0.000182027
k3 = 1.7539e-010
k4 = 5.0072346e-005
k5 = 0.000574507164
"""
Simulate this model to 200 for maximum happiness and
watch the surface to volume ratio and scaled concentrations.
"""
This example illustrates almost all of the features included
in the PySCeS MDL. Although it may be slightly more complicated
than the basic model described above it is still, by our
definition, human readable.
The PyscesUtils module holds a collection of methods of general use such as timers and
array export fucntionality that can be accessed as pysces.write.
Creates a step timer method with <name> in the TimerBox instance. Step timers
print the elapsed time as well as the next step out of maxsteps when called.
name the timer name
maxsteps the maximum number of steps associated with this timer
Constraint Based Modelling in Python (http://cbmpy.sourceforge.net)
Copyright (C) 2009-2017 Brett G. Olivier, VU University Amsterdam, Amsterdam, The Netherlands
PyscesPlot2 is a new graphics susbsystem for PySCeS which will include a
Unified Plotting Interface which can take advantage of different plotting
backends via a common user interface.
Set the size of the next plot relative to the GnuPlot canvas (e.g. screen)
size which is defined to be 1. For example if width=height=0.5 the
plot is 1/4 the size of the viewable canvas. If no arguments are supplied
reset size to 1,1.
Abstract class defining the Unified Plotting Interface methods. These methods should
be overridden and the class extended by interface specific subclasses.
This is the frontend to the PySCeS Unified Plotting Interface (pysces.plt.*) that
allows one to specify which backend should be used to plot when a
UPI method is called. More than one interface can be active at the same time and
so far the Matplotlib and GnuPlot backends are available for use.
This is an experiment which must be refactored into a more general
way of doing things. Basically, I want an instance of the abstract plotting
class which will plot to one, any or all currently available backends.
If anybody has an idea how I can generate this class automatically please let
me know ;-)
A collection of attributes defined by row and column lists
used by Response coefficients etc matrix is an array of values
while row/col are lists of row colummn name strings
Return a dictionary of <attr>_<name>, <name>_<attr> : val or {} if none
If attr exists as an index for both left and right attr then:
search=’a’ : both left and right attributes (default)
search=’l’ : left attributes only
search=’r’ : right attributes
This class is specifically designed to store the results of a time simulation
It has methods for setting the Time, Labels, Species and Rate data and
getting Time, Species and Rate (including time) arrays. However, of more use:
getOutput(*args) feed this method species/rate labels and it will return
an array of [time, sp1, r1, ….]
getDataAtTime(time) the data generated at time point “time”.
getDataInTimeInterval(time, bounds=None) more intelligent version of the above
returns an array of all data points where: time-bounds <= time <= time+bounds
getDataInTimeInterval(time, bounds=None) returns an array of all
data points where: time-bounds <= time <= time+bounds
where bound defaults to stepsize
Generic piecewise class adapted from Core2 that generates a compiled
Python code block that allows evaluation of arbitrary length piecewise
functions. Piecewise statements should be defined in assignment rules
as piecewise(<Piece>, <Conditional>, <OtherValue>) where there can
be an arbitrary number of <Piece>, <Conditional> pairs.
args a dictionary of piecewise information generated by the InfixParser as InfixParser.piecewises
Create a model object and instantiate a PySCeS model so that it can be used for
further analyses. PySCeS model descriptions can be loaded from files or strings
(see the loader argument for details).
File the name of the PySCeS input file if not explicit a *.psc extension is
assumed.
dir if specified, the path to the input file otherwise the default PyscesModel
directory (defined in the pys_config.ini file) is assumed.
autoload autoload the model, pre 0.7.1 call mod.doLoad(). (default=True) new
loader the default behaviour is to load PSC file, however, if this argument is
set to ‘string’ an input file can be supplied as the fString argument
(default=’file’)
fString a string containing a PySCeS model file (use with loader=’string’)
the File argument now sepcifies the new input file name.
Experimental: continues a simulation over a new time vector, the
CVODE memobj is reused and not reinitialised and model parameters can be
changed between calls to this method. The mod.data_sim objects from
the initial simulation and all calls to this method are stored in the list
mod.CVODE_continuous_result.
Calculate the MCA control coefficients using the current steady-state solution.
mod.__settings__[“mca_ccj_upsymb”] = 1 attach the flux control coefficients to the model instance
mod.__settings__[“mca_ccs_upsymb”] = 1 attach the concentration control coefficients to the model instance
Calculate reaction elasticities towards the parameters.
Both inputs (input1=species,input2=rates) should be valid (steady state for MCA) solutions and given in the correct order for them to be used. If either or both are missing the last state values are used automatically. Elasticities are scaled using input 1 and 2.
Arguments:
input [default=None]: species concentration vector
input2 [default=None]: reaction rate vector
Settings, set in mod.__settings__:
-elas_epar_upsymb[default=1]attachindividualelasticitysymbolstomodelinstance-elas_scaling_div0_fix[default=False]ifNaN's are detected in the variable and parameter elasticity matrix replace with zero
Calculate reaction elasticities towards the variable species.
Both inputs (input1=species,input2=rates) should be valid (steady state for MCA) solutions and given in the correct order for them to be used. If either or both are missing the last state values are used automatically.
Elasticities are scaled using input 1 and 2.
-elas_evar_upsymb[default=1]attachindividualelasticitysymbolstomodelinstance-elas_zero_conc_fix[default=False]ifzeroconcentrationsaredetectedinasteady-statesolutionmakeitaverysmallnumber-elas_zero_flux_fix[default=False]ifzerofluxesaredetectedinasteady-statesolutionmakeitaverysmallnumber-elas_scaling_div0_fix[default=False]ifINf's are detected after scaling set to zero
Forward integration steady-state solver. Finds a steady state when the
maximum change in species concentration falls within a specified tolerance.
Returns the steady-state solution and a error flag.
Algorithm controls are available as mod.fintslv_<control>
User defined forcing function either defined in the PSC input file as !F or by overwriting this method.
This method is evaluated prior to every rate equation evaluation.
PySCeS interface to the HYBRD solver. Returns a steady-state solution and
error flag. Good general purpose solver.
Algorithm controls are available as mod.hybrd_<control>
Initialise conservation related vectors/array was in InitialiseModel but has been
moved out so is can be called by when the stoichiometry is reanalysed
PySCeS interface to the LSODA integration algorithm. Given a set of initial
conditions LSODA returns an array of species concentrations and a status flag.
LSODA controls are accessible as mod.lsoda_<control>
Arguments:
initial: vector containing initial species concentrations
Initialise a PySCeS model object with PSC file that can be found in optional directory.
If a a filename is not supplied the pysces.model_dir directory contents is displayed and
the model name can be entered at the promp (<ctrl>+C exits the loading process).
Arguments:
File [default=None]: the name of the PySCeS input file
dir [default=pysces.model_dir]: the optional directory where the PSC file can be found
PySCeS interface to the (optional) NLEQ2 algorithm. This is a powerful steady-state
solver that can usually find a solution for when HYBRD() fails. Algorithm
controls are available as: mod.nleq2_<control>
Returns as steady-state solution and error flag.
PySCeS interface to the PITCON continuation algorithm. Single parameter
continuation has been implemented as a “scan” with the continuation
being initialised in mod.pitcon_par_space. The second argument does not
affect the continuation but can be used to insert a third axis parameter into
the results. Returns an array containing the results.
Algorithm controls are available as mod.pitcon_<control>
Arguments:
scanpar: the model parameter to scan (x5)
scanpar3d [default=None]: additional output parameter for 3D plots
Perform a single dimension parameter scan using the steady-state solvers. The parameter to be scanned is defined (as a model attribute “P”) in mod.scan_in while the required output is entered into the list mod.scan_out.
Results of a parameter scan can be easilly viewed with Scan1Plot().
mod.scan_in - a model attribute written as in the input file (eg. P, Vmax1 etc)
mod.scan_out - a list of required output [‘A’,’T2’, ‘ecR1_s1’, ‘ccJR1_R1’, ‘rcJR1_s1’, …]
mod.scan_res - the results of a parameter scan
mod.scan - numpy record array with the scan results (scan_in and scan_out), call as mod.scan.Vmax, mod.scan.A_ss, mod.scan.J_R1, etc.
mod.__settings__[“scan1_mca_mode”] - force the scan algorithm to evaluate the elasticities (1) and control coefficients (2)
(this should also be auto-detected by the Scan1 method).
Arguments:
range1 [default=[]]: a predefined range over which to scan.
runUF [default=0]: run (1) the user defined function mod.User_Function (!U) before evaluating the steady state.
Serialise and save a Python object using a binary pickle to file. The serialised object
is saved as <filename>.pscdat in the directory defined by mod.model_serial.
Arguments:
data: pickleable Python object
filename: the ouput filename
PySCeS integration driver routine that evolves the system over the time.
Resulting array of species concentrations is stored in the mod.data_sim object
Initial concentrations can be selected using mod.__settings__[‘mode_sim_init’]
(default=0):
0 initialise with intial concentrations
1 initialise with a very small (close to zero) value
2 initialise with results of previously calculated steady state
3 initialise with final point of previous simulation
userinit values can be (default=0):
0: initial species concentrations intitialised from (mod.S_init),
time array calculated from sim_start/sim_end/sim_points
1: intial species concentrations intitialised from (mod.S_init) existing
“mod.sim_time” used directly
2: initial species concentrations read from “mod.__inspec__”,
PySCeS non-linear solver driver routine. Solve for a steady state using HYBRD/NLEQ2/FINTSLV
algorithms. Results are stored in mod.state_species and mod.state_flux. The results
of a steady-state analysis can be viewed with the mod.showState() method.
The solver can be initialised in 3 ways using the mode_state_init switch.
mod.mode_state_init = 0 initialize with species initial values
mod.mode_state_init = 1 initialize with small values
mod.mode_state_init = 2 initialize with the final value of a 10-logstep simulation numpy.logspace(0,5,18)
Change a stoichiometric coefficient’s value in the N matrix. Only a coefficients magnitude may be set, in other words a
a coefficient’s value must remain negative, positive or zero. After changing a coefficient
it is necessary to Reanalyse the stoichiometry.
Arguments:
species: species name (s0)
reaction: reaction name (R4)
value: new coefficient value
Perform a structural analyses. The default behaviour is to construct and analyse the model
from the parsed model information. Overriding this behaviour analyses the stoichiometry
based on the current stoichiometric matrix. If load is specified PySCeS tries to load a
saved stoichiometry, otherwise the stoichiometric analysis is run. The results of
the analysis are checked for floating point error and nullspace rank consistancy.
Arguments:
override [default=0]: override stoichiometric analysis intialisation from parsed data
load [default=0]: load a presaved stoichiometry
Initialize the model stoichiometry. Given a stoichiometric matrix N, this method will return an instantiated PyscesStoich instance and status flag. Alternatively, if load is enabled, PySCeS will
attempt to load a previously saved stoichiometric analysis (saved with Stoichiometry_Save_Serial)
and test it’s correctness. The status flag indicates 0 = reanalyse stoichiometry or
1 = complete structural analysis preloaded.
Arguments:
nmatrix: The input stoichiometric matrix, N
load [default=0]: try to load a saved stoichiometry (1)
Serialize and save a Stoichiometric instance to binary pickle
Stoichiometry_Save_Serial()
Serilaise and save the current model stoichiometry to a file with name <model>_stoichiometry.pscdat
in the mod.__settings__[‘serial_dir’] directory (default: mod.model_output/pscdat)
Write an array to File with optional row/col labels. A ‘,’ separator can be specified to create
a CSV style file.
mod.__settings__[‘write_array_header’]: add <filename> as a header line (1 = yes, 0 = no)
mod.__settings__[‘write_array_spacer’]: add a space after the header line (1 = yes, 0 = no)
mod.__settings__[‘write_arr_lflush’]: set the flush rate for large file writes
Arguments:
input: the array to be written
File [default=None]: an open, writable Python file object
Row [default=None]: a list of row labels
Col [default=None]: a list of column labels
close_file [default=0]: close the file after write (1) or leave open (0)
separator [default=’ ‘]: the column separator to use
Write an array as an HTML table (no header/footer) or complete document. Tables
are formatted with coloured columns if they exceed a specified size.
mod.__settings__[‘write_array_html_header’]: write the HTML document header
mod.__settings__[‘write_array_html_footer’]: write the HTML document footer
Arguments:
input: the array to be written
File [default=None]: an open, writable Python file object
Row [default=None]: a list of row labels
Col [default=None]: a list of column labels
name [default=None]: an HTML table description line
close_file [default=0]: close the file after write (1) or leave open (0)
Write an array to an open file as a ‘LaTeX’ {array}
Arguments:
input: the array to be written
File [default=None]: an open, writable Python file object
Row [default=None]: a list of row labels
Col [default=None]: a list of column labels
close_file [default=0]: close the file after write (1) or leave open (0)
Print the moiety conserved cycles present in the system.
Arguments:
File [default=None]: an open writable Python file object
screenwrite [default=1]: write results to console (0 means no reponse)
fmt [default=’%2.3f’]: the output number format string
The PySCeS ‘save’ command, prints the entire model to screen or File in a PSC format.
(Currently this only applies to basic model attributes, ! functions are not saved).
Arguments:
filename [default=None]: the output PSC file
filepath [default=None]: the output directory
skipcheck [default=0]: skip check to see if the file exists (1) auto-averwrite
Static Bifurcation Scanning utilities using PITCON, call with loaded model object.
Hopefully nobody else was trying to use the older class as it was horrible. This new
one is is leaner, meaner and pretty cool ;-)
Performs “analysis” on the PITCON generated set of steady-state results where analysis is:
‘elasv’ = variable elasticities
‘elasp’ = parameter elasticities
‘elas’ = all elasticities
‘mca’ = control coefficients
‘resp’ = response coefficients
‘eigen’ = eigen values
‘all’ = all of the above
Higher level analysis types automatically enable the lower level analysis needed e.g.
selecting ‘mca’ implies ‘elasv’ etc. User output defined with mod.setUserOutput() is
stored in the mod.res_user array.
Arbitrary dimension generic scanner. This class is initiated with
a loaded PySCeS model and then allows the user to define scan parameters
see self.addScanParameter() and user output see self.addUserOutput().
Steady-state results are always stored in self.SteadyStateResults while
user output can be found in self.UserOutputResults - brett 2007.
While it is impossible to change the generator/range structure of a scanner
(just build another one) you can ‘in principle’ change the User Output and run
it again.
Add a parameter to scan (an axis if you like) input is:
str(name) = model parameter name
float(start) = lower bound of scan
float(end) = upper bound of scan
int(points) = number of points in scan range
bool(log) = Use a logarithmic (base10) range
bool(follower) = Scan parameters can be leaders i.e. an independent axis or a “follower” which moves synchronously with the previously defined parameter range.
Add output parameters to the scanner as a collection of one or more
string arguments (‘O1’,’O2’,’O3’, ‘On’). These are evaluated at
each iteration of the scanner and stored in the self.UserOutputResults
array. The list of output is stored in self.UserOutputList.
Returns an array of result data. I’m keepin this for backwards compatibility but
it will be replaced by a getOutput() method when this scanner is updated to use
the new data_scan object.
stst add steady-state data to output array
lbls return a tuple of (array, column_header_list)
If stst is True output has dimensions [scan_parameters]+[state_species+state_flux]+[Useroutput]
otherwise [scan_parameters]+[Useroutput].
Another looping generator function. The idea here is to create
a set of generators for the scan parameters. These generators then
all fire together and determine whether the range generators should
advance or not. Believe it or not this dynamically creates the matrix
of parameter values to be evaluated.
Evaluate the stoichiometric matrix and calculate the nullspace using LU decomposition and backsubstitution .
Generates the MCA K and Ko arrays and associated row and column vectors
Evaluate the stoichiometric matrix and calculate the left nullspace using LU factorization and backsubstitution.
Generates the MCA L, Lo, Nr and Conservation matrix and associated row and column vectors
Jordan reduction of a scaled upper triangular matrix. The returned array is now in the form [I R] and can
be used for nullspace determination. Modified row and column tracking vetors are also returned.
Core analysis algorithm; an input is preconditioned using PivotSort_initial and then cycles of PLUfactorize and
PivotSort are run until the factorization is completed. During this process the matrix is reordered by
column swaps which emulates a full pivoting LU factorization. Returns the pivot matrix P, upper factorization U
as well as the row/col tracking vectors.
Core analysis algorithm; an input is preconditioned using PivotSort_initial and then cycles of PLUfactorize and
PivotSort are run until the factorization is completed. During this process the matrix is reordered by
column swaps which emulates a full pivoting LU factorization. Returns the pivot matrix P, upper factorization U
as well as the row/col tracking vectors.
Using the R factorized form of the stoichiometric matrix we now form the K and Ko matrices. Returns
the r_ipart,Komatrix,Krow,Kcolumn,Kmatrix,Korow,info
Takes the Gauss-Jordan factorized N^T and extract the L, Lo, conservation (I -Lo) and reduced stoichiometric matrices. Returns: lmatrix_col_vector, lomatrix, lomatrix_row, lomatrix_co, nrmatrix, Nred_vector_row, Nred_vector_col, info
Arguments:
Nfull: the original stoichiometric matrix N
R_a: gauss-jordan factorized form of N^T
row_vector: row tracking vector
column_vector: column tracking vector
Performs an LU factorization using LAPACK D/ZGetrf. Now optimized for FLAPACK interface.
Returns LU - combined factorization, IP - rowswap information and info - Getrf error control.
This is a sorting routine that accepts a matrix and row/colum vectors
and then sorts them so that: there are no zero rows (by swapping with first
non-zero row) The abs(largest) pivots are moved onto the diagonal to maintain
numerical stability. Row and column swaps are recorded in the tracking vectors.
This is a sorting routine that accepts a matrix and row/colum vectors
and then sorts them so that: the abs(largest) pivots are moved onto the diagonal to maintain
numerical stability i.e. the matrix diagonal is in descending max(abs(value)).
Row and column swaps are recorded in the tracking vectors.
Calculates the dimensions of L/L0/K/K) by way of SVD and compares them to the Guass-Jordan results. Please note that for LARGE ill conditioned matrices the SVD can become numerically unstable when used for nullspace determinations
Arguments:
matrix [default=None]: the stoichiometric matrix default is self.Nmatrix
factor [default=1.0e4]: factor used to calculate the ‘zero pivot’ mask = mach_eps*factor
resultback [default=0]: return the SVD results, U, S, vh
The difference between 1.0 and the next smallest representable float
larger than 1.0. For example, for 64-bit binary floats in the IEEE-754
standard, eps=2**-52, approximately 2.22e-16.
epsnegfloat
The difference between 1.0 and the next smallest representable float
less than 1.0. For example, for 64-bit binary floats in the IEEE-754
standard, epsneg=2**-53, approximately 1.11e-16.
iexpint
The number of bits in the exponent portion of the floating point
representation.
macharMachAr
The object which calculated these parameters and holds more
detailed information.
machepint
The exponent that yields eps.
maxfloating point number of the appropriate type
The largest representable number.
maxexpint
The smallest positive power of the base (2) that causes overflow.
minfloating point number of the appropriate type
The smallest representable number, typically -max.
minexpint
The most negative power of the base (2) consistent with there
being no leading 0’s in the mantissa.
negepint
The exponent that yields epsneg.
nexpint
The number of bits in the exponent including its sign and bias.
nmantint
The number of bits in the mantissa.
precisionint
The approximate number of decimal digits to which this kind of
float is precise.
resolutionfloating point number of the appropriate type
The approximate decimal resolution of this type, i.e.,
10**-precision.
tinyfloat
The smallest positive floating point number with full precision
(see Notes).
MachAr : The implementation of the tests that produce this information.
iinfo : The equivalent for integer data types.
spacing : The distance between a value and the nearest adjacent number
nextafter : The next floating point value after x1 towards x2
For developers of NumPy: do not instantiate this at the module level.
The initial calculation of these parameters is expensive and negatively
impacts import times. These objects are cached, so calling finfo()
repeatedly inside your functions is not a problem.
Note that tiny is not actually the smallest positive representable
value in a NumPy floating point type. As in the IEEE-754 standard [1],
NumPy floating point types make use of subnormal numbers to fill the
gap between 0 and tiny. However, subnormal numbers may have
significantly reduced precision [2].
This class is specifically designed to store structural matrix information
give it an array and row/col index permutations it can generate its own
row/col labels given the label src.
Assuming that the col index array is a permutation (full/subset)
of a source label array by supplying that src to setCol
maps the row labels to cidx and creates self.col (col label list)
Assuming that the row index array is a permutation (full/subset)
of a source label array by supplying that source to setRow it
maps the row labels to ridx and creates self.row (row label list)
Set as many proxy settings as you need. You may supply a user name without
a password in which case you will be prompted to enter one (once) when required
(NO guarantees, implied or otherwise, on password security AT ALL). Arguments can be:
user = ‘daUser’,
pwd = ‘daPassword’,
host = ‘proxy.paranoid.net’,
port = 3128