Automatic Inference of Graph Transformation Rules Using the Cyclic Nature of Chemical Reactions

Abstract

Graph transformation systems have the potential to be realistic models of chemistry, provided a comprehensive collection of reaction rules can be extracted from the body of chemical knowledge. A first key step for rule learning is the computation of atom-atom mappings, i.e., the atom-wise correspondence between products and educts of all published chemical reactions. This can be phrased as a maximum common edge subgraph problem with the constraint that transition states must have cyclic structure. We describe a search tree method well suited for small edit distance and an integer linear program best suited for general instances and demonstrate that it is feasible to compute atom-atom maps at large scales using a manually curated database of biochemical reactions as an example. In this context we address the network completion problem.

Appendices

A Statistical Analysis of Rhea

Of the $M=3786$ non-isomorphic molecular graphs in RHEA, 2204 are identified uniquely by their sum formula. While 2030 of the molecules appear only in a minimum of 4 reactions, some compounds take part in a very large fraction of all reactions in RHEA, e.g., ${\mathrm{H}}^+$ participates in 11,1147 reactions, some of which are different descriptions of similar reactions where only the direction of the reaction differs, 5055 of these are truly distinct, adenosine di-, and tri-phosphate (and it derivatives), water, and dioxide each participate in more than 2000 reactions (depicted as red dots in Fig. 4 (right)). The maximum number of isomers (i.e., compounds that have the same sum formula but a non-isomorphic graph representation) is 63. The corresponding sum formula is $\mathrm{C}_15\mathrm{H}_24$. Interestingly, most of the large sets of isomers in RHEA are terpenes, condensates of identical five carbon atom building blocks. The terpenes form a combinatorial class of polycyclic ring-systems via complex sequences of cyclisation and isomerization reactions. Figure 4 (left) summarizes the results (terpenes marked with red).

Fig. 4.

Distribution of isomers and frequency of participation in reactions in Rhea. Left plot shows a few sets of isomers are very large, while most compounds in Rhea are unique up to sum formula of those compounds. Right plot shows the frequency with which a compound participates in reactions. (Color figure online)

B Analysis of Runtime

As we are mainly interested in single step reactions, we restricted our algorithms to only look for connected, vertex-disjoint transition states during the comparison. Figure 5 shows the fraction of instances where AltCyc, ILP2 and a naïve ILP-implementation with $O(n^4)$ constraints, ILP4, are able to enumerate all non-equivalent atom-atom mappings for different instance size categories as well as absolute number of instances solved divided by solution size.

Only very few instances that are not completely solved within the first 60 s are solved within reasonable time (one hour). So there seems to be a sharp divide between easy and hard instances. From the plot in Fig. 5 (left) of the fraction of instances solved fast we observe an exponential decline in ratio of solved instances. This corresponds well with the expected exponential runtime of the algorithms.

Fig. 5.

Fraction and number of instances where all optimal atom-atom maps are found in 60 s (user time) by instance size and optimal solution cost for AltCyc (magenta), ILP2 (cyan) and ILP4 (gray). (Color figure online)

As we restricted the solution set certain instances are proven infeasible by ILP2, while AltCyc will continue searching for solutions until the parameter k, the number of weight changes, gets arbitrarily high. We chose to deem instances where AltCyc found no solutions for $k\le 10$ infeasible and terminate the search. These two classes of solutions are marked in the rightmost column in Fig. 5. Note that the performance of AltCyc on the infeasible class of instances depends heavily on the somewhat arbitrary choice of maximum k.

Both ILP models are implemented using CPLEX, an efficient state of the art MIP-solver. AltCyc and ILP2 has been tested on a total of 4295 Rhea instances, while ILP4 has only been tested on a subset of these of size 250.

C Algorithmic Details

For completeness we include pseudo-code for the sub-procedures used in the paper.

Pseudo-code for WeightAlongPath: In AltCyc $^*$ (see Algorithm 2) we need to find all previous changes to an edge $\{i,j\}$ currently under examination, $w_P(\{i,j\})$.

In Algorithm 3 we show how to do this in time O(|P|), where $|P|\in O(k)$. It is possible to find $w_P(e)$ in constant time, but this would require much more complicated data structures or making changes to the graphs we work on and as k is in practice very small, this method is preferred.

To find $w_P(e)$ for a list of paths, add $w_P(e)$ for all paths in the list.

Pseudo-code for Complete : When a transition state candidate $\psi '$ is found we need to ensure it can be extended into a complete atom-atom mapping. This can be done as described in Algorithm 4. Note that the two graphs $G_1$ and $G_2$ are assumed implicitly known. The algorithm works both for a single path, P, or where P represents a list of paths.

The only non-trivial detail in Algorithm 4 is that it is not correct to remove all edges in the induced subgraph on the domain of $\psi '$, the weight change needs to be sufficient, and there may be unchanged cords to consider.

Finding 2-to-2 Candidates in $O(n^2\log n)$ Comparisons. In order to generate all $O(n^4)$ candidate reactions with no more than two molecules in the educts or products we use Algorithm 5. A set of molecules, M, is given, as well as a method to obtain the distribution of atoms and charges of the molecules h, in practice some implementation of sparse vectors. We assume we keep pointers to the original molecules that resulted in each distribution, and we get these with the function mol.

The algorithm is dominated by one of two things, either the sorting of the length $n^2$ array H (where $n = |M|$), or the time to output candidates $k\in O(n^4)$, the resulting runtime is then $O(n^2\log n + k)$.