Chemical Reaction Networks

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

Example Use

It is possible to build up a ChemicalReactionNetwork instance using the constituent parts as follows. While this is the most precise, the easiest way to create a network is using Diagrams.

First we must define some chemical species:

>>> species = sympy.var('x1 x2')
>>> species = 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, r1 represents the chemical reaction 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 ODESystem instance. The rate constants are autogenerated and follow the convention that k_i_j is the rate constant for complex i -> complex 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


A much more convenient way is to generate ChemicalReactionNetwork’s, and hence 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.x2 + 1.x1 -> 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 k_3_2 has become 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.