This project combines multiway chemical reaction networks with chemical space networks to analyze structural similarities between the products of different reaction pathways with identical or closely related precursor molecules. It explores several different avenues for analysis that can be paired with these methods, including the number of clusters in chemical space networks, the size of the largest such cluster, and branchial graph analysis. It also verifies the validity of chemical space network clustering using multiway genealogical distance and principle components analysis on the molecules produced by the multiway chemical reaction network.

Introduction

Visualization is one of the most important tools when it comes to data analysis. It allows researchers to see clear trends in the data they are working with, and different visualization techniques allow for different features of the data to be analyzed. In the past decade, chemical space networks have emerged as a powerful visualization tool for chemists and biologists. A chemical space network is a collection of molecules, each represented as a node, with an edge between two molecules if they are sufficiently similar. Chemical space networks are very versatile as they can be used on systems with a very large number of molecules, and they can be used with any similarity metric, allowing for different features of molecules to be analyzed (1).
In the world of computation, exhaustive searching is an essential tool for fully exploring a given system. In the Wolfram Language, this can be achieved with multiway graphs, which take in a list of states and iteratively apply a set of rules for a specified number of steps. The generality of multiway computation allows it to be applied in a wide array of scenarios. Previously, Ariana-Dalia Vlad completed a project developing a multiway model of chemical reactions. This allows for a great number of related molecules to be generated based on only a few reactions and a small number of starting molecules (2). In a multiway chemical system, all molecules present will be structurally related in some way, with the level of similarity depending on the reactions used. This makes the molecules from any multiway chemical reaction network a perfect candidate for analyzing using chemical space networks.
This project is heavily inspired by the work done by Robert Nachbar, especially what he presented in his 2022 Wolfram Technology Conference talk “Modeling Chemical Reaction Networks with Token-Event Graphs.” In fact, his talk presented multiway analysis on monoterpene biosynthesis, which is the main chemical system explored in this post (3).

MoleculePlotOutputForm

This is the code that automatically displays Molecule objects using MoleculePlot rather than the typical descriptor box. See the output generated below this cell for information on how to disable this temporarily.
In[]:=
moleculeImageSize[mol_Molecule,scale_:0.5]:=​​ Module[{atoms,bonds,coords,dist,minmax},​​ atoms=mol["AtomList"]/.Atom[a_,___]:>a//​​ AssociationThread[Range@Length@#,#]&;​​ bonds=Select[mol["BondList"]/.Bond[ij_,___]:>ij,​​ FreeQ[atoms/@#,"H"]&];​​ coords=QuantityMagnitude@mol["AtomDiagramCoordinates"];​​ dist=EuclideanDistance@@coords[[#]]&/@bonds//​​ Replace[{}->{0}]//Mean;​​ minmax=coords[[#]]&/@bonds//Catenate//Map[First]//MinMax;​​ ​​ Max[36,Min[6×72,(Subtract@@Reverse@minmax)/Max[1,dist]scale72]]​​ ]​​​​$RBNMoleculeFormat=True;​​​​Print["Set $RBNMoleculeFormat = False to restore built-in behavior."];​​​​Chemistry`Formatting`moleculeBox[mol_,fmt_]/;TrueQ[$RBNMoleculeFormat]:=​​ If[Chemistry`Common`has3DCoordinates[mol],​​ With[{graphics=MoleculePlot[mol]},​​ {plot=ToBoxes[graphics,fmt]},​​ Replace[plot,{HoldPattern[Graphics3DBox[x_,y___]]:>​​ Graphics3DBox[TagBox[x,BoxForm`ChemistryTag[HoldComplete[mol]]],y,​​ ContentSelectable->False,​​ Selectable->False,​​ DefaultBaseStyle->{FrontEnd`GraphicsHighlightColor->​​ RGBColor[0.269281,0.535948,0.63268]}],​​ _:>$Failed}​​ ]​​ ],​​​​With[{graphics=MoleculePlot[mol,ImageSize->moleculeImageSize[mol,0.375]]},​​ {plot=ToBoxes[graphics,fmt]},​​ Replace[plot,{HoldPattern[GraphicsBox[x_,y___]]:>​​ GraphicsBox[TagBox[x,BoxForm`ChemistryTag[HoldComplete[mol]]],y,​​ ContentSelectable->False,​​ Selectable->False,​​ DefaultBaseStyle->{FrontEnd`GraphicsHighlightColor->​​ RGBColor[0.269281,0.535948,0.63268]}],​​ _:>$Failed}​​ ]​​ ]​​ ]
Set $RBNMoleculeFormat = False to restore built-in behavior.

Pattern Reaction Paradigms

The
PatternReaction
function utilizes graph rewriting to model chemical reactions. There are several options for specifying a chemical reaction using PatternReaction, but the most powerful of these, by far, is a SMARTS string. SMARTS is an acronym that stands for SMILES (simplified molecular-input line-entry system) arbitrary target specification. The power of this method comes from allowing the user to be as specific or general as they wish to be.
For example, a condensation reaction involves the formation of a bond between two hydroxyl groups, forming an ether and giving off water as a byproduct. Utilizing SMARTS allows the user to ignore all other atoms in each molecule (or the singular molecule if it is an intramolecular condensation reaction) and only specify the atoms directly relevant to the reaction.
In[]:=
generalCondensationReaction=PatternReaction["[*:1][OH1:2].[OH1:3][*:4]>>[*:1][O:2][*:4].[OH2:3]"]
Out[]=
SMARTS also provides the ability to specify characteristics of the atoms that are not directly affected by the reaction. For example, we could specify that one of the hydroxyls is bonded to a nitrogen atom that is present in a 6-membered ring and that the other hydroxyl is bonded to a carbon atom that is present in exactly 2 rings and that these two must be in two separate molecules before the reaction takes place.
In[]:=
specificCondensationReaction=PatternReaction["([N;r6:1][OH1:2]).([OH1:3][C;R2:4])>>([N:1][O:2][C:4]).([OH2:3])"]
Out[]=
Although
PatternReaction
does not show that we have specified the ring properties for the carbon and nitrogen atoms, we can test it it on a series of molecules using the
ApplyReaction
function. Here we have 4 sets of molecules, all of which fit the pattern for
generalCondensationReaction
and only the 4th one fitting the pattern for
specificCondensationReaction
.
In[]:=
(patternReactionMoleculeSets=Map[Molecule,{{"NO","CO"},{"C1CCCCN1O","CO"},{"NO","C1CCCC2CCCCC21O"},{"C1CCCCN1O","C1CCCC2CCCCC21O"}},{2}])//Column
Out[]=

,


,


,


,

By applying
generalCondensationReaction
to all 4 sets of molecules, we expect the reactions to be carried out and receive 4 sets of reactants.
In[]:=
Row[{Plus@@#,"⟶",Plus@@ApplyReaction[generalCondensationReaction,#]},Spacer[9]]&/@patternReactionMoleculeSets//Column
Out[]=
+
⟶
+
+
⟶
+
+
⟶
+
+
⟶
+
By applying
specificCondensationReaction
to all 4 sets of molecules, we expect only the 4th set to undergo the reaction, with the first 3 reactions each returning an empty list.
In[]:=
Row[{Plus@@#,"⟶",Plus@@ApplyReaction[specificCondensationReaction,#]/.(0"∅")},Spacer[9]]&/@patternReactionMoleculeSets//Column
Out[]=
+
⟶
∅
+
⟶
∅
+
⟶
∅
+
⟶
+

Construction

Notice that we have conveniently lost track of the resultant water molecule from this condensation reaction. This allows for less processing between polymerization steps, which would help simplify code carrying out polymerization in this way.

Self-Interaction

Combination

This reaction can be represented as follows:
This method is especially useful in contexts of repeated reactions, which will be discussed in more detail in the following section on “Multiway Chemical Reaction Networks.” One of the benefits of this method is that now every reactant and product atom can be given a label and kept track of during the reaction. This helps reduce the possibility of losing or adding atoms unintentionally or putting atoms in the wrong place in the products. However, as alluded to earlier, specifying every atom in these reactions can be computationally costly for reactions involving much larger reactants.

Multiway Chemical Reaction Network

​
Note that the above function currently only works with reaction sets involving self-interaction reactions only, so more work needs to be done to generalize this to reaction sets with varying numbers of reactants.
This also provides a convenient and data-efficient format for storing all of the molecules generated by the multiway graph.

Chemical Space Network

A chemical space network is a graph where each node corresponds to an individual molecule and two nodes are connected by an edge if they have a similarity score (ranging between 0 and 1 inclusive) above a specified threshold. Depending upon the context, there are many possible choices for a similarity metric. All of the analysis in this post utilizes the TopologicalFingerprint function from WolframChemistry and compares these fingerprints using the JaccardDissimilarity function.
​
The code below is the first step in creating a chemical space network. Taking the list of SMILES strings from the multiway chemical reaction network as input, this function returns an association between pairs of molecules and their similarity scores. The pairs are stored as undirected edges to make the creation of graphs later on easier.

Statistical Analysis

There are several statistical features of the chemical space network that are of interest to graph. The three statistics I chose to analyze were:
◼
  • The number of remaining non-singleton molecules
    • ◼
  • i.e. the number of molecules similar to at least one other molecule at a given similarity threshold
    • ◼
  • The number of non-singleton connected components
    • ◼
  • The size of the largest connected component
    • The first and third statistics are monotonically decreasing as the similarity threshold is increased from 0 to 1. In general, we would like for the number of non-singletons to remain relatively high over the whole range of similarity thresholds as if too many molecules are singletons, we cannot draw as many meaningful conclusions about the similarities and differences between the molecules we are analyzing.
    Typically, the second statistic increases to a maximum at a similarity threshold slightly below 1 and then decreases a varying amount. However, there is a lot of variability in this statistic depending on the system of reactions and reactants chosen. For example, there is the possibility for several local maxima, and this seems to occur more often in smaller systems.
    Here are the plots generated for the same monoterpene system mentioned in the “Multiway Chemical Reaction Network” section, this time carried out four layers deep.
    We can see that at a similarity threshold of 1, i.e. molecules are only similar if they are topologically identical (according to their topological fingerprints), approximately one third of the molecules are singletons. The majority of the chemical reaction systems I investigated had between 10-30% of their molecules as singletons at a similarity threshold of 1.
    Here we can see that the number of connected components does not drop off as the similarity threshold approaches 1. We can also see that the chemical space network remains connected until we impose a similarity threshold slightly over 0.5 and then steadily grows to around 40 clusters.
    This graph is very similar to a step function as we can see several thresholds where the largest connected component fractionates into two or more clusters. Another behavior that can be seen with the size of the largest connected component is a steep dropoff at a high similarity threshold. based on my investigation, this seems to be more common in polymer-building systems.

    Graphs

    We would also like to display the chemical space network at various similarity thresholds to aid in our analysis of these systems. The function below creates a list of chemical space network graphs, also including the ability to click on any of the nodes to get the SMILES string of its associated molecule.
    Pictured below are the chemical space networks for the four layer deep monoterpene system described earlier at six different similarity thresholds.

    Combining the Two

    Colored Clusters on the Multiway Graph

    One way we can combine the multiway chemical reaction network and chemical space network is to display the clusters (connected components) from the chemical space network on the multiway graph using a unique color for each cluster.
    Shown below is the three layer deep monoterpene system analyzed at a similarity threshold of 0.9.
    We can see that, for the most part, molecules of the same cluster (those with the same colors) are genealogically close. There are some notable exceptions, such as light green, which has only five molecules in its cluster, but those molecules are scattered across the graph and not very genealogically close to each other.

    Genealogical Distance on Chemical Space Networks

    One other way we can combine the information from the multiway graph and the chemical space network is using the distance between the nodes in the multiway graph, which will be referred to as the genealogical distance. For each cluster in the chemical, one node is chosen as the central node, and the genealogical distances from that central node to all other nodes in the cluster are calculated. Each node is then given a color corresponding to its genealogical distance from its central node, with central nodes all colored red. The redder a node is, the closer it is genealogically to its cluster’s central node.
    Shown below is the 4 layer deep monoterpene system at a similarity threshold of 0.9 with the genealogical distance indicated by the colors.
    As we can see the image, most of the nodes within smaller clusters are genealogically very close to each other with only a few of the smaller clusters having colors besides red and orange. For the larger cluster, we see more yellow and green, which shows that the molecules within the clusters are more genealogically spread out. It is possible that these clusters would appear more red if we allowed for genealogical distance to go down trough the multiway graph, i.e. calculate the distance through “descendant” molecules, but for most analytical purposes, the definition of genealogical distance used in this analysis seems more applicable.

    Branchial Analysis

    Shown below are the branchial graphs for a two layer deep Diels-Alder reaction network discussed in the “Pattern Reaction Paradigms” section.
    Even going only two layers deep, we can see the beginnings of a cluster separation in the genealogy of the molecules. Each side of the branchial graph on the right has many more internal connections than connections with the other side.

    Full Analysis

    Example: Geranyl Cation/Monoterpene Biosynthesis Reactions (3 steps)

    Statistical analysis graphs:
    Chemical space network at a similarity threshold of 0.9:
    Multiway chemical reaction network colored using the clusters from above:
    Chemical space network with nodes colored according to genealogical distance from cluster centers:
    Principal component analysis on genealogical distance matrix from multiway graph:

    Code

    This section is a potential workflow to generate the analysis for a multiway chemical reaction system of your choice.
    The variables in this input cell are for:
    The variable in this input cell is for:
    The variable in this input cell is for:
    The variables in this input cell are for:
    This input cell returns a scatterplot of the principal components analysis on the genealogical distance matrix generated from the multiway chemical reaction network.

    Future Work

    Branchial Analysis

    As mentioned previously, branchial analysis, as it currently exists in Wolfram Language, is incompatible with most reaction networks due to the existence of directed cycles. It would be an interesting endeavor to create a version of BranchialGraph that is able to work with cyclic graphs. I suspect this would have applications far beyond the analysis done for this project, so I look forward to seeing someone tackle this problem.

    Multiway Chemical Reaction Network

    Number of Reactants

    As of now, MultiwayChemicalReactionNetwork can only handle reaction sets that have the same number of reactants for all reactions, and even then, the function has to be modified manually to accommodate different numbers of reactants. A logical next step would be to rewrite the rules given to the token event graph function in such a way as to allow for a set of reactions with different numbers of reactants.

    Pruning Function

    Many reaction systems of interest quickly explode with the number of possible molecules, slowing any analysis beyond the point of usefulness. To remedy this, a function that removes molecules at each step in the token event graph based on some input criteria would allow for larger depth investigations of some of these systems.
    There could also be a Random option for the pruning function that only allows the multiway graph to keep a certain number of molecules at each time step, choosing these at random.

    Probabilistic Multiway Computation

    One of the downsides of using multiway computation to model chemical reaction networks is that the edges in the generated graphs all have the same weights, which can be interpreted as the reactions being equiprobable. This is not very realistic, and it would be of great use to incorporate a probabilistic aspect into the multiway computation. For example, if each reaction was assigned a probability of occurring, each molecule could be paired with a list of probabilities associated with how likely it is to reach that molecule from a given state after a given number of steps. The greatest challenge with this would likely be visualization, so a great deal of creativity would be needed to incorporate this into the Wolfram Language.

    Chemical Space Network

    Other Similarity Metrics

    The analysis with the current implementation can only be done using the topological fingerprints and the JaccardDissimilarity function for similarity comparisons. To expand the analytical abilities of this method, a more general form for the ChemicalSpaceNetwork function allowing for different metrics to be used should be implemented

    GUI

    A WSS ’23 project by Joshua Khorsandi explored chemical space networks in the context of comparing bio activity to chemical structure. As part of his project, he developed a GUI that allowed for individual cluster investigation and molecular inspection among other capabilities (4). His work should be modelled after and expanded upon to incorporate the multiway chemical reaction network and analyses described in this post.

    Concluding Remarks

    The versatility of multiway computation has great potential in chemical analysis, especially when paired with more traditional analysis methods, such as chemical space networks. From the results presented, it is clear that multiway chemical reaction networks and chemical space networks provide very similar, yet distinct, information about the chemical systems that they are analyzing. This makes pairing them together a very powerful tool: verifying the validity of each other while also providing more information than either method could on its own.

    Acknowledgments

    I would like to acknowledge Robert Nachbar, who served as my mentor throughout this project. He consistently helped me focus my project towards goals that were both significant and achievable in the time frame of the Wolfram Summer School and swiftly resolve any coding challenges I ran into during my project. I would also like to thank Stephen Wolfram and Mads Bahrami for directed me towards this project and providing their insights.

    References

    1
    .
    G. Maggiora and J. Bajorath (2014), “Chemical space networks: a powerful new paradigm for the description of chemical space,” J Comput Aided Mol Des, https://pubmed.ncbi.nlm.nih.gov/24925682/
    2
    .
    [FSWS23] A Multiway Model of Chemical Reactions, by Ariana-Dalia Vlad, Wolfram Community, January 13, 2023, https://community.wolfram.com/groups/-/m/t/2778062
    3
    .
    Modeling Chemical Reaction Networks with Token-Event Graphs, by Robert Nachbar, Wolfram Technology Conference, October 21, 2022, https://www.wolfram.com/broadcast/video.php?c=104&p=6&disp=list&v=3812
    4
    .
    Chemical space networks from molecule data, by Joshua Khorsandi, Wolfram Community, July 12, 2023, https://community.wolfram.com/groups/-/m/t/2958091

    Cite This Notebook

    “Analyzing multiway chemical reaction networks using chemical space networks”
    by Nicholas Frieler​
    ​https://community.wolfram.com/groups/-/m/t/3210073