Fast Enumeration of Non-isomorphic Chemical Reaction Networks

Abstract

Chemical reaction networks (CRNs) have been applied successfully to model a wide range of phenomena and are commonly used for designing molecular computation circuits. Often, CRNs with specific properties (oscillations, Turing patterns, multistability) are sought, which entails searching an exponentially large space of CRNs for those that satisfy a property. As the size of the CRNs being considered grows, so does the frequency of isomorphisms, by up to a factor $N!$, where $N$ is the number of species. Accordingly, being able to generate sets of non-isomorphic CRNs within a class can lead to large computational savings when carrying out global searches. Here, we present a bijective encoding of bimolecular CRNs into novel vertex-coloured digraphs called . The problem of enumerating non-isomorphic CRNs can then be tackled by leveraging well-established computational methods from graph theory [20]. In particular, we extend Nauty, the graph isomorphism tool suite by McKay [22]. Our method is highly parallelisable and more efficient than competing approaches, and a software package (genCRN) is freely available for reuse. Non-isomorphs are generated directly by genCRN, alleviating the need to store intermediate results. We provide the first complete count of all 2-species bimolecular CRNs and extend previous known counts for classes of CRNs of special interest, such as mass-conserving and reversible CRNs.

Appendices

A Definitions

This section introduces definitions for automorphisms, orbits and the autormophism group for graphs, following [31].

Definition 6 (Permutation)

A permutation of a set S is a total function from S to itself.

Definition 7 (Cyclic permutation)

A permutation $\pi$ of the form:

$$ \left( \begin{array}{l l l l l l} x &{} \pi (x) &{} \pi ^2(x) &{} \cdots &{} \pi ^{p-2}(x) &{} \pi ^{p-1}(x) \\ \pi (x) &{} \pi ^2(x) &{} \pi ^3(x) &{} \cdots &{} \pi ^{p-1}(x) &{} x \end{array} \right) $$

is said to be cyclic permutation of period p.

Definition 8 (Disjoint cycle representation)

A disjoint cycle representation of a permutation $\pi$ on a set S is a composition of cyclic permutations on subsets of S that constitute a partition of S, one cyclic permutation for each subset in the partition.

Definition 9 (Group)

An algebraic system ${<}U, \star {>}$ is a called a group if it has the following properties:

• 1. the operation $\star$ is associative,

• 2. there is an identity element,

• 3. every element of U has an inverse.

Definition 10 (Permutation group)

A closed non-empty collection P of permutations on a set Y of objects that forms a group under the operation of composition is called a permutation group. The combined structure may be denoted $V = [P : V]$. It is often denoted P when the set of Y objects is understood from context.

Definition 11 (Orbit)

Let $\mathcal {P} = [P:Y]$ be a permutation group, and let $y \in Y$. The orbit of the object y under the action of P is the set $\{\pi (y) ~|~ \pi \in P \}$.

Corollary 1

Let $\mathcal {P} = [P:Y]$ be a permutation group. Then being coorbital is an equivalence relation

Proof

Identity: by the identity permutation. Commutativity: because each $\pi$ is invertible. Transitivity: by function composition $\circ$.

B Proofs

Lemma 1

Let $\mathcal {N}$ be a bimolecular CRN. Then $\llbracket \mathcal {N}\rrbracket ^{I_1}_{J_1} \cong \llbracket \mathcal {N}\rrbracket ^{I_2}_{J_2}$ holds for any indexing sets $I_1, I_2, J_1, J_2$.

Proof

The lemma can be proved by explicitly constructing bijections $\alpha$ and $\beta$ required by Definition 5. Recall that we indicate with $\{c_{i}\}_{i \in I}$ and $\{A_{j}\}_{j \in J}$ respectively the indexed set of the complexes and of the species in $\mathcal {N}$.

Let $\alpha = \{(i_1, i_2)~|~ c_{i_1} = c_{i_2} \text { for } i_1 \in I_1, i_2 \in I_2\}$ and $\beta = \{(j_1, j_2)~|~ A_{j_1} = A_{j_2} \text { for } j_1 \in J_1, j_2 \in J_2 \}$. These functions are well-defined because the indexing sets all target the same CRN $\mathcal {N}$. It is also easy to show that they are bijections.

The lemma is proved by verifying that $\alpha$ and $\beta$ satisfy Condition 1 to 3 of Definition 5:

• 1. $V_1\alpha \beta = V_2$ because $\alpha$ and $\beta$ are bijections over the indexed sets;

• 2. $E_1\alpha \beta = \{(j_1, i_1)~|~A_{j_1} \!\in \! c_{i_1}\}\alpha \beta \,\cup \, \{(i_1, i_1')~|~ c_{i_1} \rightarrow c_{i_1'} \in \mathcal {R}\}\alpha \beta$   by Definition 4

  $= \{(j_2, i_2)~|~A_{j_2} \in c_{i_2}\} \cup \{(i_2, i_2')~|~ c_{i_2} \rightarrow c_{i_2'} \in \mathcal {R}\}$    by def. of $\alpha$, $\beta$.

  $= E_2$

  which proves the case.

• 3. Let $i_1$ be such that $\sigma _1(i_1) = \emptyset$. By Definition 4, $c_{i_1} = \emptyset$, and since $\alpha (i_1) = i_2$ implies that $c_{i_1} = c_{i_2}$, then $c_2 = \emptyset$ as well. Therefore $\sigma (i_2) = \emptyset$ holds by Definition 4, which implies $\sigma (i_1)\alpha = \sigma _2(i_2)$. The proof for the remaining cases (monomers, homodimers and heterodimers) is similar.    $\square$

Lemma 2

Let $\mathcal {N}_1$ and $\mathcal {N}_2$ be bimolecular CRNs. If $\mathcal {N}_1 \cong \mathcal {N}_2$, then $\llbracket \mathcal {N}_1\rrbracket ^{I_1}_{J_1} \cong \llbracket \mathcal {N}_2\rrbracket ^{I_2}_{J_2}$ for any indexing sets $I_1, I_2, J_1$ and $J_2$.

Proof

By definition of CRN isomorphism (Definition 2), there exists a permutation $\pi$ over the species of $\mathcal {N}_1$ such that $\mathcal {N}_1\pi = \mathcal {N}_2$. Notice that by Lemma 1 we can deduce $\llbracket \mathcal {N}_1\pi \rrbracket ^{I_1}_{J_1} \cong \llbracket \mathcal {N}_2\rrbracket ^{I_2}_{J_2}$ for any indexing sets $I_1, I_2, J_1$ and $J_2$. The proof of this lemma is similar to Lemma 1, by defining $\alpha = \{(i_1, i_2)~|~ c_{i_1}\pi = c_{i_2} \text { for } i_1 \in I_1, i_2 \in I_2\}$ and $\beta = \{(j_1, j_2)~|~ A_{j_1}\pi = A_{j_2} \text { for } j_1 \in J_1, j_2 \in J_2 \}$.    $\square$

Lemma 3

Let $\mathcal {N}_1$ and $\mathcal {N}_2$ be bimolecular CRNs with indexing sets respectively $I_1, J_1$ and $I_2, J_2$. If $\llbracket \mathcal {N}_1\rrbracket ^{I_1}_{J_1} \cong \llbracket \mathcal {N}_2\rrbracket ^{I_2}_{J_2}$, then $\mathcal {N}_1 \cong \mathcal {N}_2$.

Proof

Let $\llbracket \mathcal {N}_1\rrbracket ^{I_1}_{J_1} = \langle V_1, E_1, \sigma _1, \rho _1\rangle$ and $\llbracket \mathcal {N}_2\rrbracket ^{I_2}_{J_2} = \langle V_2, E_2, \sigma _2, \rho _2 \rangle$, such that $\llbracket \mathcal {N}_1\rrbracket ^{I_1}_{J_1} \cong \llbracket \mathcal {N}_2\rrbracket ^{I_2}_{J_2}$. By hypothesis, $\llbracket \mathcal {N}_1\rrbracket ^{I_1}_{J_1} \cong \llbracket \mathcal {N}_2\rrbracket ^{I_2}_{J_2}$ implies the existence of bijections $\alpha$ and $\beta$ that satisfy conditions 1 to 3 in Definition 5. Let us define the following permutation of $\mathcal {S}$:

$$\begin{aligned} \pi&= \{(A_{j_1}, A_{j_2})~|~ j_1\beta = j_2 \} \circ \pi _I \end{aligned}$$

where $\pi _I$ is the identity function over $\mathcal {S}$. Since $\beta$ and $\pi _I$ are bijections, then $\pi$ is also a bijection; since its domain and range are $\mathcal {S}$, $\pi$ is a well-defined permutation.

Let $c_{i_1} \rightarrow c_{i_1'}$ be a reaction in $\mathcal {N}_1$. By Definition 4 $E_1$ contains the edge $(i_1, i_1')$. Since $\llbracket \mathcal {N}_1\rrbracket ^{I_1}_{J_1}$ and $\llbracket \mathcal {N}_1\rrbracket ^{I_2}_{J_2}$ are isomorphic by hypothesis, it follows by definition that $E_2 = E_1\alpha \beta$, therefore the edge $(i_1, i_1')\alpha = (i_2, i_2')$ also exists in $E_2$ for $i_2, i_2' \in I_2$. Because of this, the reaction $c_{i_2} \rightarrow c_{i_2'}$ exists in $\mathcal {N}_2$; moreover, by Condition 3 of Definition 5, the complexes have the same stoichiometry.

Similarly, let $(j_1, i_1)$ be an edge in $E_1$ such that $A_{j_1}$ occurs in $c_{i_1}$. By Condition 2 of Definition 5, $E_2$ contains the edge $(j_1, i_1)\alpha \beta = (j_1\alpha , i_1\beta ) = (i_2, j_2)$, which means that $A_{j_2}$ occurs in $c_{i_2}$. By definition of $\pi$, $A_{j_1}\pi = A_{j_1\beta } = A_{j_2}$; since $c_{i_1}$ and $c_{i_2}$ also have the same multiplicity by Condition 3 of Definition 5, this implies that $c_{i_1}\pi = c_{i_2}$. A similar line of reasoning shows that $c_{i_1'}\pi = c_{i_2'}$. Therefore $(c_{i_1} \rightarrow c_{i_1'})\pi = c_{i_2} \rightarrow c_{i_2'}$. Generalising this result to all reactions in $\mathcal {N}_1$, we obtain $\mathcal {N}_1 \pi = \mathcal {N}_2$, which concludes the proof.    $\square$

Lemma 4

The number of p-CRNs (reactions have up to p reactants/products) with up to $N$ species grows as $\mathcal {O}(2^{N^{2p}})$.

Proof

Following Eq. 3, the total number of p-CRNs with up to $N$ and specifically $M$ reactions is given by ${ L_p(N)(L_p(N)-1) \atopwithdelims ()M}$. Given that

$$\begin{aligned} L(N)(L(N)-1) = \dfrac{(N+p)\dots (N+1)}{p!} . \dfrac{(N+p)\dots (N+1) - p!}{p!} = \mathcal {O}(N^{2p}), \end{aligned}$$

we can use the fact that $\sum _{i=1}^k{n\atopwithdelims ()k}=2^n$ to characterise the total number of bimolecular CRNs as $\mathcal {O}(2^{N^{2p}})$.    $\square$

C Complex-Multiplicity-Species Graph

Section 2.1 has shown how to encode bimolecular CRNs as vertex-coloured digraphs. It is natural to wonder whether this encoding extends to more than bimolecularity. Unfortunately cannot encode higher molecularities than 2, however we propose in this section a more general encoding of CRNs called the Complex-Multiplicity-Species graph (or CMS graph).

We begin by showing that cannot encode trimolecular reactions. Consider in fact the reaction $2A + B \rightarrow A + 2B$. If we added a new color “$2\square +\square$” and connect two node species A and B to it, there would be no way to tell which of the two species is actually the homodimer and which one is the monomer.

Fig. 9.

Complex-Multiplicity-Species graph  encoding of a CRN.

To overcome this issue, we propose Complex-Multiplicity-Species graphs, which extend with multiplicity nodes, that is distinct coloured nodes that point out the multiplicity of a species in a reaction. If m is the molecularity of interest, then there are $m+1$ kinds of multiplicity nodes: naught, $\square$, $2\square$, $3\square$ and so on. Each species node is connected to m multiplicity nodes, signifying for example A, 2A, 3A etc. Naught is a separate multiplicity node that cannot be connected to any species. In turn, multiplicity nodes are connected to complex nodes to represent the original CRN’s reaction. Figure 9 show an example of a CMS graph; notice that no confusion is possible between complexes $2A + B$ and $A + 2B$.

We believe that CMS graphs are a general encoding of CRNs with any molecularity, but we leave a formal definition and proofs for future work.

D Counts of Non-isomorphic CRNs

In this appendix, we tabulate the numbers of non-isomorphic CRNs found using genCRN. The tables can be compared against values reported at https://reaction-networks.net/networks/, at the time of writing, which were evaluated using the method in [3]. In each case, we report values for genuine CRNs, those which use all $N$ species.

1.1 D.1 No Filters

Here, we consider the total number of non-isomorphic CRNs for $N$ species and $M$ reactions. The results are graphically depicted in Fig. 5, but tabulated below (Table 1).

Table 1. Genuine non-isomorphic CRNs. The number of non-isomorphic CRNs is shown for different numbers of species and reactions. Coloured in blue are those counts not available at https://reaction-networks.net/networks/ at time of writing.

1.2 D.2 Reversible CRNs

To generate reversible CRNs, we generate undirected graphs of a suitable size and feed these into genCRNin the same way as for general CRNs with irreversible reactions. Reported below are counts for $M$ reversible reactions. As such, the CRNs found have $2M$ unidirectional reactions.

Fig. 10.

Counts and execution times for enumeration of genuine non-isomorphic reversible CRNs.

Table 2. Genuine non-isomorphic reversible CRNs. The number of non-isomorphic CRNs with only reversible reactions is shown for different numbers of species and reactions. Coloured entries correspond to comparisons with the counts available at https://reaction-networks.net/networks/ at time of writing. Blue indicates values not available, and red indicates values that differ.

1.3 D.3 Non-trivial Dynamics

As described in the main text, one can test whether a CRN has non-trivial dynamics. To apply this filter to the enumerated non-isomorphic CRNs, one can use the -t flag for genCRN.

Fig. 11.

Counts and execution times for enumeration of genuine non-isomorphic CRNs with non-trivial dynamics.

Table 3. Genuine non-isomorphic CRNs with non-trivial dynamics. The number of non-isomorphic CRNs with non-trivial dynamics is shown for different numbers of species and reactions. Coloured entries correspond to comparisons with the counts available at https://reaction-networks.net/networks/ at time of writing. Blue indicates values not available, and red indicates values that differ.

We found that our counts differ with those reported at https://reaction-networks.net/networks/ (using the method in [3]) for CRNs with at least 4 species and 4 reactions (Table 3). In each of the 5 counts identified as different, the values we report are lower than those reported previously, though within 1.5% relative error. By analysing the sets of CRNs produced in the 6 species and 4 reactions class, we have determined that we incorrectly classify at least one CRN as trivial, and tracked this to numerical discrepancies in our use of Fourier-Motzkin elimination. Further work would be required to refine the numerical procedure used here, in order to improve confidence in our counts of dynamically non-trivial and non-isomorphic CRNs. For instance, implementing Fourier-Motzkin elimination with explicit rational numbers (one integer variable each for the numerator and denominator), could avoid the emergence of values close to zero.