The pybrn Tutorial

Package initialization:

>>> import brn

It is useful to also load the numpy package:

>>> import numpy as np

The brn package provides the following routines:

brn.Brn

Constructor for a new biochemical network object. The Brn class provides the main interface for functionality of the pybrn package, via the class’ methods.

brn.fromSBML

Import a network object from an SBML file.

brn.fromSBMLString

Import a network object from an SBML string. Useful for example for directly downloading network definitions from the internet, see doc/examples/biomodels_direct.py for an example script.

brn.make_massaction_network

Generate a network with irreversible mass action reaction rates from a given stoichiometric matrix.

First, a few other modules are required for the purpose of this tutorial:

>>> import os, brn.test.generic

Creation of a biochemical network object:

Create a network directly from code:

>>> net = brn.Brn([[1,-1]],species=['x'],reactions=['v1','v2'],parameters=['k'], rates={'v1':2.5, 'v2':'k*x'},parvals={'k':0.1})

Import from SBML:

>>> netsbml = brn.fromSBML(os.path.join("doc","examples","simplenet.xml"))

Generate from stoichiometric matrix (assuming irreversible mass action reactions):

>>> netma = brn.make_massaction_network([[2,-2],[1,-1]])

Handling model data

Relevant model data, may be given as numbers or symbolic expressions (strings, python syntax): - parameters - initial conditions for species concentrations - reaction rate expressions - variables

Setting model data with Brn.set, for example:

>>> net.set({'k':0.5, 'x':1, 'v1':'2*k'})

All numerical values and symbolic expressions for model data are stored in the dict ``Brn.values´´, and can be accessed there:

>>> parameter_value = net.values['k']
>>> print parameter_value
0.5

Model evaluation

The Brn class supports evaluation of arbitrary symbolic expressions in the context of a given Brn object via the method Brn.eval. The method can also be used to obtain numerical values for the network’s initial conditions or parameter values, for example:

>>> p0 = net.eval(net.parameters)
>>> x0 = net.eval(net.species)

This also offers a convenient way to evaluate for example reaction rates. The valuedict optional argument can be used to give custom numerical values for species concentrations or parameters:

>>> v0 = net.eval(net.reactions,valuedict={'my_species':10,'my_parameter':0.1})

Conservation relations

Checking whether a network has conserved moieties:

>>> cons_net = brn.test.generic.make_conservation_net()
>>> print cons_net.has_conservation
True

Do an analysis of conservation relations via SVD

The method Brn.conservation_analysis computes a reduced stoichiometric matrix, link matrix, and vectors describing the conservation relations. These are stored in the attribute Brn.conserv:

• Brn.conserv.Nr: reduced stoichiometric matrix
• Brn.conserv.L: link matrix
• Brn.conserv.G: matrix of conservation relations

Example for such an analysis:

>>> conservation_data = cons_net.conservation_analysis()
>>> np.dot(cons_net.conserv.G.T,cons_net.stoich)        # should be equal to zero
array([[  1.11022302e-16,   1.11022302e-16]])
>>> np.dot(cons_net.conserv.L,cons_net.conserv.Nr)      # shoud be equal to net.stoich
array([[ -1.00000000e+00,  -8.89045781e-17],
       [  1.00000000e+00,  -1.00000000e+00],
       [  1.04083409e-17,   1.00000000e+00]])

functions to map from reduced network with variable z to original network with variable x:

>>> x0 = [3,0,0]
>>> z2x = cons_net.conserv.reduced2full_fun(x0)
>>> x2z = cons_net.conserv.full2reduced_fun()
>>> z2x(x2z(x0))
array([  3.00000000e+00,  -1.11022302e-16,   2.22044605e-16])

Steady state computation

Steady states are computed via the Brn.steadystate method:

>>> xs = net.steadystate()

It is also possible to give custom parameter values:

>>> xs = net.steadystate(values={'k':0.5})

The computation is a local search which is initalized at the default species concentration values for the network. In addition, a custom starting point can be given as optional parameter in the method call:

>>> xs = net.steadystate(x0=[1])

If you want to compute the steady states for a range of parameter values, use the Brn.steadystate_branch method. The following call computes the steady states under variation of the parameter k from between values 0.5 and 1.5:

>>> plist, xslist = net.steadystate_branch('k',0.5,1.5)

As in the steadystate method, custom parameter values can be given via the optional argument values:

>>> plist, xslist = net.steadystate_branch('k',0.5,1.5,values={'v1':5})

Simulation

The Brn.simulate method computes a trajectory of the network’s dynamics, returning state values at time points as given in the method’s argument:

>>> tspan = np.arange(101.)/10
>>> t,x = net.simulate(tspan)
Found integrator vode

As with steady state computations, custom parameter values and initial conditions can be given in the values optional argument:

>>> tc,xc = net.simulate(tspan,values={'x':0.1,'k':0.2})
Found integrator vode

A custom initial condition, for example from a previous simulation, can also be given in the optional argument x0:

>>> t,x = net.simulate(np.arange(tspan[-1],2*tspan[-1],0.1),x0=xc[-1])
Found integrator vode