Chemical Reaction Networks
==========================

One flavour of dynamical system that can be analysed using desr is the Chemical Reaction Network (CRN).
The :py:class:`~desr.chemical_reaction_network.ChemicalReactionNetwork` class is useful for drawing out dynamical systems and automatically generating :py:class:`~desr.ode_system.ODESystem` instances.

Example Use
-----------

It is possible to build up a :py:class:`~desr.chemical_reaction_network.ChemicalReactionNetwork` instance using the constituent parts as follows.
While this is the most precise, the easiest way to create a network is using :ref:`diagram_example`.

First we must define some chemical species:

    >>> species = sympy.var('x1 x2')
    >>> species = list(map(ChemicalSpecies, species))
    >>> x1, x2 = species

Then we must make complexes, which are representations of one side of a chemical reaction diagram.
It is represented using a dictionary, where the keys represent the chemical elements and the value represents its coefficient in the equation.

    >>> complex0 = Complex({x1: 1, x2: 1})
    >>> complex1 = Complex({x2: 2})
    >>> complex2 = Complex({x1: 1})
    >>> complex3 = Complex({x2: 1})
    >>> complexes = (complex0, complex1, complex2, complex3)


The final step is to make reactions, which simply take two complexes (two sides of an equation) and form an arrow from the first side to the second.
For example, :code:`r1` represents the chemical reaction :code:`x1 + x2 -> 2*x2`.

    >>> r1 = Reaction(complex0, complex1)
    >>> r2 = Reaction(complex3, complex2)
    >>> reactions = [r1, r2]

Finally, we form a network from these constituent parts and form an :py:class:`~desr.ode_system.ODESystem` instance.
The rate constants are autogenerated and follow the convention that :code:`k_i_j` is the rate constant for complex :code:`i` -> complex :code:`j`.

    >>> reaction_network = ChemicalReactionNetwork(species, complexes, reactions)
    >>> system = reaction_network.to_ode_system()
    >>> system
    dt/dt = 1
    dx1/dt = -k_0_1*x1*x2 + k_3_2*x2
    dx2/dt = k_0_1*x1*x2 - k_3_2*x2
    dk_0_1/dt = 0
    dk_3_2/dt = 0
    >>> answer = ODESystem.from_equations('dx1/dt = -k_0_1*x1*x2 + k_3_2*x2\ndx2/dt = k_0_1*x1*x2 - k_3_2*x2')
    >>> system == answer
    True

.. _diagram_example:

Diagrams
--------

A much more convenient way is to generate :class:`~desr.chemical_reaction_network.ChemicalReactionNetwork`'s, and hence :class:`~desr.ode_system.ODESystem`'s, with diagrams.

This can be done like so:

    >>> reactions = ['x1 + x2 -> 2*x2',
    ...              'x2 -> x1']
    >>> reaction_network = ChemicalReactionNetwork.from_diagram('\n'.join(reactions))
    >>> reaction_network
    1.x1 + 1.x2 -> 2.x2
    1.x2 -> 1.x1
    >>> system = reaction_network.to_ode_system()
    >>> system
    dt/dt = 1
    dx1/dt = -k_0_1*x1*x2 + k_2_3*x2
    dx2/dt = k_0_1*x1*x2 - k_2_3*x2
    dk_0_1/dt = 0
    dk_2_3/dt = 0

Note that :code:`k_3_2` has become :code:`k_2_3` since complex order is automatically determined when feeding in a diagram.
Otherwise, we have an entirely equivalent way of constructing CRNs and their associated dynamical systems.