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