sign_crn.reaction_networks

Define species, complexes, and reaction networks.

(Chemical) reaction networks

This module provides tools for defining species, complexes, and (chemical) reaction networks with (generalized) mass-action kinetics. It includes functionality for analyzing reaction networks, such as checking weak reversibility, computing deficiencies, and verifying conditions for existence and uniqueness of positive complex-balanced equilibria (CBE) based on [MHR19] and [AMR24].

We define species for a reaction network:

sage: from sign_crn import *
sage: A, B, C = species("A, B, C")

Next, we create a reaction network and add complexes and reactions:

sage: rn = ReactionNetwork()
sage: rn.add_complex(0, A + B)
sage: rn.add_complex(1, C)
sage: rn.add_reactions([(0, 1), (1, 0)])
sage: rn.plot()
Graphics object consisting of 6 graphics primitives

We analyze the reaction network:

sage: rn.is_weakly_reversible()
True
sage: rn.deficiency_stoichiometric
0
sage: rn.has_exactly_one_cbe()
True

For further examples, see the class ReactionNetwork.

Functions

species(names)

Define species from a string of names.

Classes

Complex(species_dict)	A complex involving species.
ReactionNetwork()	A reaction network with (generalized) mass-action kinetics.

class sign_crn.reaction_networks.Complex(species_dict: Dict[_Species, int | float | var])

A complex involving species.

This class represents a linear combination of species, supporting various operations such as addition, subtraction, and scalar multiplication. It also supports symbolic expressions and substitution of values for variables.

EXAMPLES:

First, we define some species:

sage: from sign_crn import *
sage: species("A, B, C")
(A, B, C)

Usual operations like addition and multiplication are supported:

sage: A + B
A + B
sage: 2 * A + 3 * B
2*A + 3*B

Symbolic expressions work as well:

sage: var("a")
a
sage: 2 * a * A
2*a*A
sage: (2 + a) * A
(a + 2)*A
sage: a * (A + B)
a*A + a*B

Similar to polynomials, we can substitute values for variables in a complex:

sage: complex = a * A - B
sage: complex
a*A - B
sage: complex(a=0)
-B
sage: complex(a=1)
A - B
sage: (a * A)(a=0)
0
sage: A(a=1)
A

TESTS:

Operations with invalid types raise appropriate errors:

sage: A * B
Traceback (most recent call last):
...
TypeError: Cannot multiply species by species.
sage: B + A
A + B
sage: A * 2
2*A
sage: 0 * A
0
sage: 1 * A
A
sage: A + 0
A
sage: A + 1
Traceback (most recent call last):
...
TypeError: Cannot add <class 'sage.rings.integer.Integer'> to species.
sage: -A
-A
sage: (-a - 1) * A
(-a - 1)*A
sage: A - B
A - B
sage: A - 2 * B
A - 2*B
sage: (2 * A + 3 * B).get_coefficient(A)
2

LaTeX representation is supported for better visualization:

sage: species("A, B")
(A, B)
sage: 2 * A + 3 * B
2*A + 3*B
sage: (2 * A + 3 * B)._latex_()
'2 \\, A + 3 \\, B'
sage: (A - B)._latex_()
'A - B'
sage: (A + B)._latex_()
'A + B'
sage: (2 * A - 3 * B)._latex_()
'2 \\, A - 3 \\, B'
sage: (A)._latex_()
'A'
sage: (0 * A)._latex_()
'0'

static from_species(species: _Species) → Complex

Return a complex from a species.

get_coefficient(species: _Species | Complex) → int | float

Return the coefficient of the species in the complex.

If a species is not present, return 0.

involved_species() → set[_Species]

Return the species involved in the complex.

class sign_crn.reaction_networks.ReactionNetwork

A reaction network with (generalized) mass-action kinetics.

This class represents a reaction network, where complexes are connected by directed reactions. It supports generalized mass-action kinetics and symbolic rate constants.

This class provides tools for:

• Adding and removing complexes and reactions.

• Computing stoichiometric and kinetic-order matrices.

• Analyzing network properties, such as weak reversibility and deficiencies.

• Checking conditions for uniqueness and existence of positive complex-balanced equilibria (CBE).

• Visualizing the reaction network as a directed graph.

EXAMPLES:

We define a reaction network with two complexes involving variables in the kinetic orders:

sage: from sign_crn import *
sage: var("a, b")
(a, b)
sage: species("A, B, C")
(A, B, C)
sage: rn = ReactionNetwork()
sage: rn.add_complex(0, A + B, a * A + b * B)
sage: rn.add_complex(1, C)
sage: rn.add_reactions([(0, 1), (1, 0)])
sage: rn
Reaction network with 2 complexes, 2 reactions and 3 species.
sage: rn.complexes_stoichiometric
{0: A + B, 1: C}
sage: rn.complexes_kinetic_order
{0: a*A + b*B, 1: C}
sage: rn.reactions
[(0, 1), (1, 0)]
sage: rn.species
(A, B, C)
sage: rn.plot()
Graphics object consisting of 6 graphics primitives

We describe the reaction network using matrices:

sage: rn.matrix_of_complexes_stoichiometric
[1 0]
[1 0]
[0 1]
sage: rn.matrix_of_complexes_kinetic_order
[a 0]
[b 0]
[0 1]

The stoichiometric and kinetic-order matrices are given by:

sage: rn.stoichiometric_matrix
[-1  1]
[-1  1]
[ 1 -1]
sage: rn.kinetic_order_matrix
[-a  a]
[-b  b]
[ 1 -1]

We check some conditions for our reaction network:

sage: rn.are_both_deficiencies_zero()
True
sage: rn.is_weakly_reversible()
True
sage: rn(a=2, b=1).has_robust_cbe()
True
sage: rn.has_robust_cbe()
[{a > 0, b > 0}]
sage: rn.has_at_most_one_cbe()
[{a >= 0, b >= 0}]

To extend our network, we add further complexes and reactions:

sage: var("c")
c
sage: species("D, E")
(D, E)
sage: rn.add_complexes([(2, D, c * A + D), (3, A), (4, E)])
sage: rn.add_reactions([(1, 2), (3, 4), (4, 3)])
sage: rn
Reaction network with 5 complexes, 5 reactions and 5 species.
sage: rn.plot()
Graphics object consisting of 15 graphics primitives

To make this network weakly reversible, we add another reaction:

sage: rn.is_weakly_reversible()
False
sage: rn.add_reaction(2, 0)
sage: rn.is_weakly_reversible()
True

Now, our network consists of 6 reactions:

sage: rn.reactions
[(0, 1), (1, 0), (1, 2), (2, 0), (3, 4), (4, 3)]

The corresponding rate constants are:

sage: rn.rate_constants()
(k_0_1, k_1_0, k_1_2, k_2_0, k_3_4, k_4_3)

We compute the incidence and source matrices of the corresponding directed graph:

sage: rn.incidence_matrix()
[-1  1  0  1  0  0]
[ 1 -1 -1  0  0  0]
[ 0  0  1 -1  0  0]
[ 0  0  0  0 -1  1]
[ 0  0  0  0  1 -1]
sage: rn.source_matrix()
[1 0 0 0 0 0]
[0 1 1 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 0]
[0 0 0 0 0 1]

The Laplacian matrix involving the rate constants is given by:

sage: rn.laplacian_matrix()
[        -k_0_1          k_1_0          k_2_0              0              0]
[         k_0_1 -k_1_0 - k_1_2              0              0              0]
[             0          k_1_2         -k_2_0              0              0]
[             0              0              0         -k_3_4          k_4_3]
[             0              0              0          k_3_4         -k_4_3]
sage: rn.ode_rhs()
(-k_0_1*x_A^a*x_B^b + k_2_0*x_A^c*x_D - k_3_4*x_A + k_1_0*x_C + k_4_3*x_E,
 -k_0_1*x_A^a*x_B^b + k_2_0*x_A^c*x_D + k_1_0*x_C,
 k_0_1*x_A^a*x_B^b - (k_1_0 + k_1_2)*x_C,
 -k_2_0*x_A^c*x_D + k_1_2*x_C,
 k_3_4*x_A - k_4_3*x_E)

The network is described by the following matrices:

sage: rn.matrix_of_complexes_stoichiometric
[1 0 0 1 0]
[1 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
[0 0 0 0 1]
sage: rn.matrix_of_complexes_kinetic_order
[a 0 c 1 0]
[b 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
[0 0 0 0 1]
sage: rn.stoichiometric_matrix
[-1  1  0  1 -1  1]
[-1  1  0  1  0  0]
[ 1 -1 -1  0  0  0]
[ 0  0  1 -1  0  0]
[ 0  0  0  0  1 -1]
sage: rn.kinetic_order_matrix
[   -a     a     c a - c    -1     1]
[   -b     b     0     b     0     0]
[    1    -1    -1     0     0     0]
[    0     0     1    -1     0     0]
[    0     0     0     0     1    -1]
sage: rn.stoichiometric_matrix_as_kernel
[1 0 1 1 1]
[0 1 1 1 0]
sage: rn.kinetic_order_matrix_as_kernel
[    1     0     a a - c     1]
[    0     1     b     b     0]

We check some conditions for our network:

sage: rn.deficiency_stoichiometric
0
sage: rn.deficiency_kinetic_order
0
sage: rn.is_weakly_reversible()
True
sage: rn(a=2, b=1, c=1).has_robust_cbe()
True
sage: rn.has_robust_cbe() # random order
[{a > 0, a - c > 0, b > 0}]
sage: rn(a=2, b=1, c=1).has_at_most_one_cbe()
True
sage: rn.has_at_most_one_cbe() # random order
[{a >= 0, a - c >= 0, b >= 0}]
sage: rn.has_exactly_one_cbe()
Traceback (most recent call last):
...
ValueError: Method does not support parameters in the complexes.
sage: rn(a=2, b=1, c=1).has_exactly_one_cbe()
True

Next, we remove one component and a reaction of our network:

sage: rn.remove_complex(3)
sage: rn.remove_complex(4)
sage: rn.remove_reaction(1, 0)
sage: rn
Reaction network with 3 complexes, 3 reactions and 4 species.

We define a reaction network involving molecules:

sage: A, B, C = species("H_2, O_2, H_2O")
sage: var("a")
a
sage: rn = ReactionNetwork()
sage: rn.add_complex(0, 2 * A + B, 2 * a * A + a * B)
sage: rn.add_complex(1, 2 * C)
sage: rn.species
(H_2, H_2O, O_2)
sage: rn.add_reactions([(0, 1), (1, 0)])
sage: rn.plot()
Graphics object consisting of 6 graphics primitives

We can also define an ecosystem involving animals as species:

sage: fox, rabbit = species("Fox, Rabbit")
sage: rn = ReactionNetwork()
sage: rn.add_complex(0, rabbit + fox)
sage: rn.add_complex(1, 2 * fox)
sage: rn.add_reactions([(0, 1), (1, 0)])
sage: rn.plot()
Graphics object consisting of 6 graphics primitives

add_complex(i: int, complex_stoichiometric: Complex, complex_kinetic_order: Complex | None = None) → None

Add complex to reaction network.

add_complexes(complexes: List[Tuple[int, Complex, Complex | None]]) → None

Add complexes to reaction network.

add_reaction(start: int, end: int) → None

Add reaction to reaction network.

add_reactions(reactions: List[Tuple[int, int]]) → None

Add reactions to reaction network.

are_both_deficiencies_zero() → bool

Return whether both deficiencies are zero.

property deficiency_kinetic_order: int

Return the kinetic-order deficiency.

property deficiency_stoichiometric: int

Return the stoichiometric deficiency.

has_at_most_one_cbe() → bool

Check whether there is at most one positive CBE in every stoichiometric class and for all rate constants.

has_exactly_one_cbe() → bool

Check whether there is a unique positive CBE in every stoichiometric class and for all rate constants.

Note

This method does not support symbolic expressions (variables) in the complexes.

has_exactly_one_equilibrium() → bool

Check whether there is a unique positive equilibrium.

Note

This method does not support symbolic expressions (variables) in the complexes.

EXAMPLES:

sage: from sign_crn import *
sage: species("X, Y")
(X, Y)
sage: rn = ReactionNetwork()
sage: rn.add_complexes([(1, 2 * X + Y), (2, 2 * Y), (3, 3 * X + Y), (4, 4 * X), (5, 3 * X + 2 * Y), (6, 3 * X + 3 * Y), (7, 2 * X + 4 * Y)])
sage: rn.add_reactions([(1, 2), (3, 4), (5, 6), (6, 7)])
sage: rn.has_exactly_one_equilibrium()
True
sage: rn.remove_reaction(6, 7)
sage: rn.has_exactly_one_equilibrium()
True
sage: rn.remove_reaction(1, 2)
sage: rn.has_exactly_one_equilibrium()
False

sage: rn = ReactionNetwork()
sage: rn.add_complexes([(0, 5 * X + Y), (1, 6 * X), (2, 4 * X + 3 * Y), (3, 5 * X + 6 * Y), (4, 2 * X + 3 * Y), (5, 4 * Y), (6, 4 * X + 2 * Y), (8, 3 * X + 3 * Y), (9, 2 * X + 5 * Y)])
sage: rn.add_reactions([(0, 1), (2, 3), (4, 5), (6, 2), (8, 9)])
sage: rn.has_exactly_one_equilibrium()
True
sage: rn.remove_reaction(0, 1)
sage: rn.has_exactly_one_equilibrium()
False

has_robust_cbe() → bool

Check whether there is a unique positive CBE in every stoichiometric class, for all rate constants and for all small perturbations of the kinetic orders.

incidence_matrix(**kwargs) → matrix

Return the incidence matrix of the graph.

is_weakly_reversible() → bool

Return whether each component of the reaction network is strongly connected.

property kinetic_order_matrix: matrix

Return the kinetic-order matrix where the columns correspond to the reactions.

Each columns stands for a reaction, and each row stands for a species.

property kinetic_order_matrix_as_kernel: matrix

Return the kernel matrix of the kinetic-order matrix.

laplacian_matrix() → matrix

Return the Laplacian matrix of the graph.

property matrix_of_complexes_kinetic_order: matrix

Return the matrix that decodes the kinetic-order complexes of the reaction network.

Each column stands for a complex, and each row stands for a species.

property matrix_of_complexes_stoichiometric: matrix

Return the matrix that decodes the stoichiometric complexes of the reaction network.

Each column stands for a complex, and each row stands for a species.

ode_rhs() → vector

Return the right hand side of the ordinary differential equation of this reaction network.

plot(kinetic_orders: bool = True, layout: str = 'spring', edge_labels: bool = False, vertex_colors: str | List[str] = 'white', vertex_size: int = 6000, **kwargs) → None

Plot the reaction network.

This method visualizes the reaction network as a directed graph. The vertices represent complexes, and the edges represent reactions. The layout, labels, and other visual properties can be customized.

INPUT:

• kinetic_order (bool, default: True): If True, displays both stoichiometric and kinetic-order complexes for each vertex. If False, only stoichiometric complexes are shown.

• layout (str, default: “spring”): Specifies the layout of the graph. Common options include “circular” and “spring”.

• edge_labels (bool, default: False): If True, displays the rate constants as labels on the edges.

• vertex_colors (str or list, default: “white”): Specifies the color of the vertices. Can be a single color or a list of colors corresponding to each vertex.

• vertex_size (int, default: 6000): Specifies the size of the vertices in the plot.

• **kwargs: Additional keyword arguments passed to the underlying graph plotting function.

OUTPUT:

• A graphical representation of the reaction network.

EXAMPLES:

sage: from sign_crn import *
sage: species("A, B, C")
(A, B, C)
sage: rn = ReactionNetwork()
sage: rn.add_complex(0, A + B)
sage: rn.add_complex(1, C)
sage: rn.add_reactions([(0, 1), (1, 0)])
sage: rn.plot()
Graphics object consisting of 6 graphics primitives

We can customize plotting:

sage: rn.plot(edge_labels=True)
Graphics object consisting of 8 graphics primitives
sage: rn.plot(kinetic_orders=False)
Graphics object consisting of 6 graphics primitives
sage: rn.plot(vertex_size=3000)
Graphics object consisting of 6 graphics primitives
sage: rn.plot(layout="circular")
Graphics object consisting of 6 graphics primitives

rate_constants() → Tuple[var, ...]

Return rate constants.

property reactions: List[Tuple[int, int]]

Return reactions.

remove_complex(i: int) → None

Remove complex from reaction network.

remove_reaction(start: int, end: int) → None

Remove reaction from reaction network.

set_rate_constant_variable(variable: var) → None

Set rate constant variable. This method allows you to set a custom variable for the rate constants.

EXAMPLES:

sage: from sign_crn import *
sage: species("A, B, C")
(A, B, C)
sage: rn = ReactionNetwork()
sage: rn.add_complexes([(0, A + B), (1, C)])
sage: rn.add_reactions([(0, 1), (1, 0)])
sage: rn.rate_constants()
(k_0_1, k_1_0)
sage: rn.plot(edge_labels=True)
Graphics object consisting of 8 graphics primitives

You can also use a variable with a LaTeX name:

sage: rn.set_rate_constant_variable(var("tau"))
sage: rn.rate_constants()
(tau_0_1, tau_1_0)
sage: rn.plot(edge_labels=True)
Graphics object consisting of 8 graphics primitives
sage: var("k", latex_name=r"\kappa")
k
sage: rn.set_rate_constant_variable(k)
sage: rn.rate_constants()
(k_0_1, k_1_0)
sage: rn.plot(edge_labels=True)
Graphics object consisting of 8 graphics primitives

source_matrix(**kwargs) → matrix

Return the source matrix of the graph.

property species: Tuple[Complex, ...]

Return the species of the reaction network as a tuple of complexes.

property stoichiometric_matrix: matrix

Return the stoichiometric matrix where the columns correspond to the reactions.

Each columns stands for a reaction, and each row stands for a species.

property stoichiometric_matrix_as_kernel: matrix

Return the kernel matrix of the stoichiometric matrix.

sign_crn.reaction_networks.species(names: str) → Complex | Tuple[Complex, ...]

Define species from a string of names.

The string can contain species names separated by commas or spaces. The function defines each species globally, similar to var().

See Complex for operations and more details.

EXAMPLES:

sage: from sign_crn import *
sage: species("A")
A
sage: A
A
sage: species("A, B, C")
(A, B, C)
sage: A
A
sage: B
B
sage: C
C
sage: species("A,B,C")
(A, B, C)
sage: species("A B C")
(A, B, C)
sage: A, B = species("H_2, O_2")
sage: A
H_2
sage: B
O_2