sign_vector_conditions.reaction_networks

Setting up reaction networks

This module provides tools for defining species, complexes, and 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 unique positive complex-balanced equilibria (CBE).

Key Classes and Functions:

• ReactionNetwork: Represents a reaction network.

• species(): Utility function for defining species.

For detailed examples and usage, see ReactionNetwork.

EXAMPLES:

We define species for a reaction network:

sage: from sign_vector_conditions import *
sage: A, B, C = species("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 analyze the reaction network:

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

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_vector_conditions.reaction_networks.Complex(species_dict: Dict[sign_vector_conditions.reaction_networks._Species, Union[int, float, sage.calculus.var.var]])

A complex involving species.

This class represents a linear combination of species with coefficients, 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_vector_conditions 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 are also supported:

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: sign_vector_conditions.reaction_networks._Species) → sign_vector_conditions.reaction_networks.Complex

Return a complex from a species.

get_coefficient(species: Union[sign_vector_conditions.reaction_networks._Species, sign_vector_conditions.reaction_networks.Complex]) → Union[int, float]

Return the coefficient of the species in the complex.

If a species is not present, return 0.

involved_species() → set[sign_vector_conditions.reaction_networks._Species]

Return the species involved in the complex.

class sign_vector_conditions.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, allowing for symbolic rate constants and kinetic orders.

The ReactionNetwork 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 unique positive complex-balanced equilibria (CBE).

• Visualizing the reaction network as a directed graph.

Key Attributes:

• graph: The directed graph representing the reaction network.

• complexes_stoichiometric: A dictionary mapping complex indices to stoichiometric complexes.

• complexes_kinetic_order: A dictionary mapping complex indices to kinetic-order complexes.

• species: A tuple of all species involved in the network.

Key Methods:

• add_complex(), remove_complex(): Add or remove complexes from the network.

• add_reaction(), remove_reaction(): Add or remove reactions between complexes.

• plot(): Visualize the reaction network as a directed graph.

• stoichiometric_matrix(), kinetic_order_matrix(): Get the stoichiometric and kinetic-order matrices.

• has_robust_cbe(), has_at_most_one_cbe(), has_exactly_one_cbe(): Check conditions for unique positive CBE.

EXAMPLES:

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

sage: from sign_vector_conditions 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]
sage: rn.stoichiometric_matrix_as_kernel
[1 0 1]
[0 1 1]
sage: rn.kinetic_order_matrix_as_kernel
[1 0 a]
[0 1 b]

We check some conditions for our system:

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}]

We extend our network by adding 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 system 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 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.differential_equation()
(-k_0_1*x_0^a*x_2^c*x_3 + k_1_0*x_0^b + k_2_0*x_1 - k_3_4*x_2 + k_4_3*x_4,
 -k_0_1*x_0^a*x_2^c*x_3 + k_1_0*x_0^b + k_2_0*x_1,
 k_0_1*x_0^a*x_2^c*x_3 - (k_1_0 + k_1_2)*x_0^b,
 k_1_2*x_0^b - k_2_0*x_1,
 k_3_4*x_2 - k_4_3*x_4)

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 system:

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 variables!
sage: rn(a=2, b=1, c=1).has_exactly_one_cbe()
True

We remove one component and a reaction of our system:

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.

Here is an example 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: sign_vector_conditions.reaction_networks.Complex, complex_kinetic_order: Optional[sign_vector_conditions.reaction_networks.Complex] = None) → None

Add complex to system.

add_complexes(complexes: List[Tuple[int, sign_vector_conditions.reaction_networks.Complex, Optional[sign_vector_conditions.reaction_networks.Complex]]]) → None

Add complexes to system.

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

Add reaction to system.

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

Add reactions to system.

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.

differential_equation() → sage.modules.free_module_element.vector

Return the differential equation of the system.

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_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) → sage.matrix.constructor.matrix

Return the incidence matrix of the graph.

is_weakly_reversible() → bool

Return whether each component of the system is strongly connected.

property kinetic_order_matrix: sage.matrix.constructor.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: sage.matrix.constructor.matrix

Return the kernel matrix of the kinetic-order matrix.

laplacian_matrix() → sage.matrix.constructor.matrix

Return the Laplacian matrix of the graph.

property matrix_of_complexes_kinetic_order: sage.matrix.constructor.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: sage.matrix.constructor.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.

plot(kinetic_orders: bool = True, layout: str = 'spring', edge_labels: bool = False, vertex_colors: Union[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_vector_conditions 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[sage.calculus.var.var, ...]

Return rate constants.

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

Return reactions.

remove_complex(i: int) → None

Remove complex from system.

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

Remove reaction from system.

set_rate_constant_variable(variable: sage.calculus.var.var) → None

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

EXAMPLES:

sage: from sign_vector_conditions 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) → sage.matrix.constructor.matrix

Return the source matrix of the graph.

property species: Tuple[sign_vector_conditions.reaction_networks.Complex, ...]

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

property stoichiometric_matrix: sage.matrix.constructor.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: sage.matrix.constructor.matrix

Return the kernel matrix of the stoichiometric matrix.

sign_vector_conditions.reaction_networks.species(names: str) → Union[sign_vector_conditions.reaction_networks.Complex, Tuple[sign_vector_conditions.reaction_networks.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 as var().

See Complex for operations and more details.

EXAMPLES:

sage: from sign_vector_conditions 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