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
True
Diagrams¶
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.