Modeling Chemical Reaction Networks with Token-Event Graphs

Robert B. Nachbar, Wolfram Solutions

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

The recent addition of PatternReaction and ApplyReaction to Wolfram Language now allows the user to programmatically manipulate the atoms and bonds in Molecule objects in a chemically meaningful way. This talk will show how these functions can be used with token-event graphs to model chemical reaction networks. Examples with the Diels-Alder reaction, Cope rearrangement, Wagner-Meerwein rearrangement, and terpene biosynthesis will be demonstrated. As these networks represent multiway systems, they can also be studied from that perspective with branchial and causal graphs.

Outline

Introduction

Reaction Networks

Metabolic Networks

Organic Synthesis

Rearrangements

Models of Chemistry

Molecules

Reactions

Models for Network Generation

Multiway Systems

Chemistry with TEG—Vinyl Chloride Polymerization

Caveats

Chain Reaction

Diels-Alder Reaction

Rearrangements

Cope Rearrangement

Wagner-Meerwin Rearrangement

Carbocation Rearrangements

Norbornyl Cation Rearrangement

Metabolism

Terpene Biosynthesis

Strategy

More PatternReactions

Monoterpene Biosynthesis Token-Event Graph

Conclusions

Acknowledgements

Initialization

Introduction

Reaction Networks

Metabolic Networks

One of the most common examples of a reaction network is the biochemistry of metabolism

Boehringer Ingelheim (now Roche) first published this 4' × 6' wall chart in the 1970’s, and it’s now available on the web

Module[{detailedView,schematicOverview,tcaCycle,lipidMetabolism, cite},
detailedView=;
schematicOverview=;
tcaCycle=;
lipidMetabolism=;
cite="http://biochemical-pathways.com/#/map/1";
TabView[{"Full Details"->citation[Show[detailedView,ImageSize->Large],cite],
"Schematic Overview"->citation[Show[schematicOverview,ImageSize->Large],cite],
"TCA Cycle"->citation[Show[tcaCycle,ImageSize->Large],cite],
"Lipd Metabolism"->citation[Show[lipidMetabolism,ImageSize->Large],cite]},
ControlPlacement->Left
]
]
centerAndCloseCellGroup[]

Organic Synthesis

Another common example is the multiple step route to complex organic molecules

Many syntheses are linear, beginning with the starting material and incrementally adding pieces (usually carbon atoms) until the product is attained

Woodward's synthesis of strychnine

The total synthesis of strychnine was first published by Robert B. Woodward in 1954, and there are at least 11 subsequently published syntheses

Module[{cite},
cite="https://en.wikipedia.org/wiki/Strychnine_total_synthesis";
TabView[{
"1 to 9"->citation[Show[Import["https://upload.wikimedia.org/wikipedia/commons/thumb/4/45/Strychnine_Woodward_1954_start.svg/1112px-Strychnine_Woodward_1954_start.svg.png"],ImageSize->600],cite],
"9 to 15"->citation[Show[Import["https://upload.wikimedia.org/wikipedia/commons/thumb/0/04/Strychnine_Woodward_part_2.svg/1028px-Strychnine_Woodward_part_2.svg.png"],ImageSize->600],cite],
"15 to 20"->citation[Show[Import["https://upload.wikimedia.org/wikipedia/commons/thumb/8/85/Strychnine_Woodward_part_3.svg/1012px-Strychnine_Woodward_part_3.svg.png"],ImageSize->600],cite],
"20 to 23"->citation[Show[Import["https://upload.wikimedia.org/wikipedia/commons/thumb/2/27/Strychnine_Woodward_part_4.svg/956px-Strychnine_Woodward_part_4.svg.png"],ImageSize->600],cite],
"23 to 28"->citation[Show[Import["https://upload.wikimedia.org/wikipedia/commons/thumb/a/a9/Strychnine_Woodward_part_5.svg/928px-Strychnine_Woodward_part_5.svg.png"],ImageSize->600],cite]
},
ControlPlacement->Left
]]
centerAndCloseCellGroup[]

Others are convergent, first building two or more large portions of the final product and then joining them together near the end

Corey's synthesis of aspidophytine

Module[{cite},
cite="https://synarchive.com/syn/38";
TabView[{
"1 to 6"->citation[Show[,ImageSize->16 72],cite],
"7 to 15"->citation[Show[,ImageSize->16 72],cite],
"6 + 15 to 22"->citation[Show[,ImageSize->16 72],cite]
},
ControlPlacement->Left
]]
centerAndCloseCellGroup[]

Rearrangements

Some organic molecules are metastable and under certain conditions will rearrange the bonds among the atoms, whose number remains constant

The Cope rearrangement is an electrocyclic reaction, changing the bonds via a cyclic intermediate

citation[Show[Import["https://upload.wikimedia.org/wikipedia/commons/c/c8/Cope.png"] // ImageTake[#,80]&,ImageSize->600],"https://en.wikipedia.org/wiki/Cope_rearrangement"]
centerAndCloseCellGroup[]

The Wagner-Meerwein rearrangement 1,2-shifts of hydrogen or carbon atoms

citation[Show[Import["https://upload.wikimedia.org/wikipedia/commons/thumb/4/46/Wagner-Meerwein_Camphen_mechanism.svg/2000px-Wagner-Meerwein_Camphen_mechanism.svg.png"],ImageSize->16 72],"https://en.wikipedia.org/wiki/Wagner–Meerwein_rearrangement"]
centerAndCloseCellGroup[]

Models of Chemistry

So, how does one get from these nice drawings to something computable?

There are two aspects of chemistry that we need to consider:

Structure: the atoms and how they are connected to each other via chemical bonds (pairs of electrons)

Reactions: changing the connections

We will ignore energetics (energy barriers of transition states) and mechanistic details (single step vs. multiple steps)

Let’s start with an example: reduction of citronellal with sodium borohydride to citronellol:

Models are used to represent both molecules and their reactions

George E. P. Box: “All models are wrong, but some are useful”

Under the hood, the molecular structure is treated as a mathematical graph (where the colored vertices represent the atoms and the labeled edges represent the bonds), and the reaction is treated as graph rewriting (where the edges are modified, added, or deleted)

The numbers of kinds of atoms and the number of electron pairs are conserved

Molecules

Molecules are are represented in WL, naturally enough, by the symbol Molecule

Here is the reactant molecule, citronellal:

\[KernelIcon];(citro=Molecule["citronellal"])//

Reactions

The reaction involves only a few of the atoms of citronellal, so we can abstract away most of the molecule with an “R-group”

TabView[{
"Reactants"->,
Column[{"arrow pushing ="," graph rewriting"}]->With[{img=},ImageResize[ImageTake[img,{36,All}],720]],
"Products"->
},ControlPlacement->Left,FrameMargins->36]
centerAndCloseCellGroup[]

The graph rewriting rules are made with the symbol PatternReaction

We’ll make use of the third form that uses SMARTS notation to define the reaction:

\[KernelIcon];rxn=PatternReaction["([C:1][C:2](=[O:3])[H:4]).([Na+:10].[H:5][B-:6]([H:7])([H:8])[H:9]).([H:11][O:12][C:13])>>([C:1][C:2]([O:3][H:11])([H:4])[H:5]).([B+0:6]([H:7])([H:8])[H:9]).([Na+:10].[O-:12][C:13])",ImageSize->1200]

The graph rewriting is carried out with the function ApplyReaction

These are the reagents that participate in the reaction:

\[KernelIcon];sodiumBorohidride=Molecule["sodium borohydride"]

\[KernelIcon];methanol=Molecule["methanol"]

Finally, the reaction can be evaluated with

\[KernelIcon];products=ApplyReaction[rxn,{citro,sodiumBorohidride,methanol}]

\[KernelIcon];MapThread[MoleculePlot[#,ImageSize->#2]&,{products,{Medium,Small,Small}}]//Row[#,Style[" + ",72]]&
centerAndCloseCellGroup[]

In this example, the reaction had only one place in the starting material at which it could be applied, so it was “done” after this one “step”

In reaction networks there are many opportunities for reaction, and in the next section we’ll describe one way they can be explored

Models for Network Generation

Multiway Systems

A multiway system consists of an initial state, a set of rules that can modify the state, and a procedure that applies the rules in all possible ways to a state to generate new states

It’s an iterative process

It may or may not terminate

There are many kinds of multiway systems; one of the simplest is string rewriting

We can use ResourceFunction["MultiwaySystem"]

\[KernelIcon];ResourceFunction["MultiwaySystem"][{"A"->"AB","B"->"A"},"A",5,"EvolutionGraph"]//LayeredGraphPlot

Or we can use ResourceFunction["TokenEventGraph"], which is more general

TEG uses the resource function, and is defined in the initialization section at the end of the notebook

\[KernelIcon];TEG[s_:>StringReplaceList[s,#]&/@{"A"->"AB","B"->"A"},"A",5]

Chemistry with TEG—Vinyl Chloride Polymerization

Vinyl chloride polymerization is a chain reaction that can be studied with token-event graphs

The first few steps are shown here

Module[{peroxidePattern,alkoxyRadicalPattern,initiation,radicalPattern,vinylPattern,coupledRadicalPattern,elongation,peroxide,vinylChloride,products1,products2,products3,products4,plus,arrow,miniscule},
peroxidePattern=MoleculePattern[{Atom["C","DisplayName"->"R"],Atom["O"],Atom["O"],Atom["C","DisplayName"->"R"]},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Single"]}];
alkoxyRadicalPattern=MoleculePattern[{Atom["C","DisplayName"->"R"],Atom["O","UnpairedElectronCount"->1,"DisplayName"->"(O^•)"]},{Bond[{1,2},"Single"]}];

initiation=PatternReaction[{peroxidePattern}->{alkoxyRadicalPattern,alkoxyRadicalPattern}](*//Echo*);

radicalPattern=MoleculePattern[{Atom["C"|"O","UnpairedElectronCount"->1,"DisplayName"->"(O|C)^•"]},{}];
vinylPattern=MoleculePattern[{"C","C","Cl"},{Bond[{1,2},"Double"],Bond[{2,3}]}];
coupledRadicalPattern=MoleculePattern[{Atom[_(*,"DisplayName"->"(O|C)"*)],Atom["C"],Atom["C","UnpairedElectronCount"->1,"DisplayName"->"(C^•)"],Atom["Cl"]},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Single"]}];

elongation=PatternReaction[{radicalPattern,vinylPattern}->{coupledRadicalPattern},{{1,1}->{1,1},{2,1}->{1,2},{2,2}->{1,3},{2,3}->{1,4}}](*//Echo*);

peroxide=Molecule[{Atom["C","DisplayName"->"R"],Atom["O"],Atom["O"],Atom["C","DisplayName"->"R"]},{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Single"]}];
vinylChloride=Molecule["C=CCl"];

products1=ApplyReaction[initiation,{peroxide}];
products2=ApplyReaction[elongation,{First@products1,vinylChloride}];
products3=ApplyReaction[elongation,{First@products2,vinylChloride}];
products4=ApplyReaction[elongation,{First@products3,vinylChloride}];

plus=Style["+",36];
arrow=Style["⟶",36];
miniscule=45;
TabView[{
"radical formation"->Row[{MoleculePlot[peroxide,ImageSize->Tiny],Row[MoleculePlot[#,ImageSize->miniscule]&/@products1,plus]},arrow],
"initiation"->Row[{Row[{MoleculePlot[First@products1,ImageSize->miniscule],MoleculePlot[vinylChloride,ImageSize->Tiny]},plus],Row[MoleculePlot[#,ImageSize->Small]&/@products2,plus]},arrow],
"elongation"->Row[{Row[{MoleculePlot[First@products2,ImageSize->Small],MoleculePlot[vinylChloride,ImageSize->Tiny]},plus],Row[MoleculePlot[#,ImageSize->Small]&/@products3,plus]},arrow],
"elongation"->Row[{Row[{MoleculePlot[First@products3,ImageSize->Small],MoleculePlot[vinylChloride,ImageSize->Tiny]},plus],Row[MoleculePlot[#,ImageSize->Medium]&/@products4,plus]},arrow]
},
ControlPlacement->Left,
FrameMargins->36
]
]
centerAndCloseCellGroup[]

This code defines the elongation reaction of the polymerization process:

\[KernelIcon];radicalPattern=MoleculePattern[];
vinylPattern=MoleculePattern[];
coupledRadicalPattern=MoleculePattern[];
elongation=PatternReaction[{radicalPattern,vinylPattern}->{coupledRadicalPattern},{{1,1}->{1,1},{2,1}->{1,2},{2,2}->{1,3},{2,3}->{1,4}}]

Here are the initial alkoxide radical and vinyl chloride:

\[KernelIcon];alkoxideRadical=Molecule[];
vinylChloride=Molecule[];

We can test the polymerization with NestList:

\[KernelIcon];NestList[ApplyReaction[elongation,{Last@#,vinylChloride}]&,{alkoxideRadical},5]//Catenate

And now we can generate the token-event graph model for vinyl chloride polymerization:

\[KernelIcon];TEG[radical_:>ApplyReaction[elongation,{radical,vinylChloride}],alkoxideRadical,5,"TokenRenderingFunction"->(MoleculePlot[#,]&),]

The token vertices (molecules) are rendered with MoleculePlot, and the event vertices (reactions) are rendered as filled yellow circles

Caveats

It should be noted that multiway systems are a mathematical construct, and our use here is not meant to represent a simulation of a real physical system

In general, we will ignore quantitative aspects of geometry, stereochemistry, and energetics

On the other hand, we also don’t want to get so caught up in mathematical rigor that we loose sight of the chemistry

Chain Reaction

Diels-Alder Reaction

The Diels-Alder reaction joins a conjugated diene molecule with an alkene or alkyne molecule in an electrocyclic condensation

m1=Molecule["C=CC=C",AtomDiagramCoordinates->];
m2=Molecule["C=CC",AtomDiagramCoordinates->];
m3=Molecule["CC1CC=CCC1",AtomDiagramCoordinates->];
Grid[{{MoleculePlot[m1,IncludeHydrogens->False,ImageSize->54],Style[" +  ",24],MoleculePlot[m2,IncludeHydrogens->False,ImageSize->54],Style[" ⟶ ",24],MoleculePlot[m3,IncludeHydrogens->False,ImageSize->144]},{"diene","","dienophile","","product"}},BaseStyle->{"Text",24,ShowStringCharacters->False}]
centerAndCloseCellGroup[]

The pattern reaction can be written conveniently with SMARTS:

\[KernelIcon];dielsAlder=PatternReaction["[C:1]=[C:2][C:3]=[C:4].[C:5]=[C:6]>>[C:1]1[C:2]=[C:3][C:4][C:5][C:6]1"]

We can demonstrate this reaction with butadiene and propene

\[KernelIcon];butadiene=Molecule[];
propene=Molecule[];

\[KernelIcon];ApplyReaction[dielsAlder,{butadiene,propene}]//

Notice that the product contains a double bond, and therefore can act as a dienophile, thus setting the stage for a chain reaction

The structural diversity of the products becomes more interesting if we use these reactants:

\[KernelIcon];{cyclohexadiene=Molecule["C1=CC=CCC1"],
isoprene=Molecule["C=C(C)C=C"]}//

Each of these molecules can act as either the diene or the dienophile!

\[KernelIcon];ApplyReaction[dielsAlder,{cyclohexadiene,isoprene},All]//Catenate//

\[KernelIcon];ApplyReaction[dielsAlder,{isoprene,cyclohexadiene},All]//Catenate//

Here is one step of the token-event graph

\[KernelIcon];TEG[{{diene_,dienophile_}:>Catenate@ApplyReaction[dielsAlder,{diene,dienophile},All]},{cyclohexadiene,isoprene},1,"TokenDeduplication"->True,"TokenLabeling"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenRenderingFunction"->,ImageSize->20 72]

Examination of the products reveals that none of them contain conjugated diene moieties, so the process will terminate after one step

Requesting two steps demonstrates this condition:

\[KernelIcon];TEG[{{diene_,dienophile_}:>Catenate@ApplyReaction[dielsAlder,{diene,dienophile},All]},{cyclohexadiene,isoprene},2,"TokenDeduplication"->True,"TokenLabeling"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenRenderingFunction"->,ImageSize->20 72]

However, we can supply the necessary conjugated dienes with the "TokenDecompositionFunction"->(DeleteDuplicates[Join[#,{cyclohexadiene,isoprene}],MoleculeMatchQ]&) option

\[KernelIcon];g=TEG[{{diene_,dienophile_}:>Catenate@ApplyReaction[dielsAlder,{diene,dienophile},All]},
{},2,
"TokenDecompositionFunction"->(DeleteDuplicates[Join[#,{cyclohexadiene,isoprene}],MoleculeMatchQ]&),
"EventDeduplication"->True,
"TokenDeduplication"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,
"TokenLabeling"->True,"TokenRenderingFunction"->,
ImageSize->30 72,AspectRatio->1/4
]
closeCellGroup[]

The molecules highlighted in yellow are present in the PubChem or the Chemical Abstracts Service Registry databases

The database membership has been determined on the fly with MoleculeValue

\[KernelIcon];MoleculeValue[isoprene,"PubChemCompoundID"]

\[KernelIcon];MoleculeValue[isoprene,"CASRegistryNumber"]

It’s somewhat surprising that none of the molecules from the second step are known

Here are some stats for this token-event graph:

{{"# Reactions",Count[VertexList[g],{"Event",__}]},
{"# Molecules",Count[VertexList[g],_Molecule]},
{"# Duplicates",VertexInDegree[g,#]&/@Cases[VertexList[g],_Molecule]//Count[#,2]&}}//
Grid[#,Alignment->{{Left,Right}}]&
closeCellGroup[]

These are the molecules that are generated by more than one route

\[KernelIcon];duplicates=Cases[VertexList[g],m_Molecule/;VertexInDegree[g,m]==2]

These structures can be grouped by similarity: {{1,2,3,4},{5,6},{7},{8,9,10,11},{12}}

Here are the “synthetic routes” for molecules 7 and 12

\[KernelIcon];extractRouteDAR[{Molecule[],Molecule[]}]//
centerAndCloseCellGroup[]

Here are the “synthetic routes” for molecules 8, 9, 10, and 11

\[KernelIcon];extractRouteDAR[{Molecule[],Molecule[],Molecule[],Molecule[]},ImageSize->4 72]//
centerAndCloseCellGroup[]

Seeding the system with cyclohexadiene and allowing only isoprene for subsequent condensations (as studied by Andersen at al.) gives

\[KernelIcon];TEG[{{mol_}:>Catenate@Join[ApplyReaction[dielsAlder,{isoprene,mol},All],ApplyReaction[dielsAlder,{mol,isoprene},All]]},{cyclohexadiene},4,"EventVertices"->False,"TokenDeduplication"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenLabeling"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenRenderingFunction"->,GraphLayout->"RadialEmbedding",ImageSize->18 72,AspectRatio->1]
centerAndCloseCellGroup[]

The event (reaction) vertices have been omitted for clarity

Rearrangements

Cope Rearrangement

The Cope rearrangement is a [3,3]-sigmatropic rearrangement reaction of a 1,5-diene

m1=Molecule["CC=CCCC=C",AtomDiagramCoordinates->];
m2=MoleculeModify[m1,{"DeleteBond",{{2,3},{3,4},{4,5},{5,6},{6,7}}}]//MoleculeModify[#,{"AddBond",{Bond[{2,3},"Single"],Bond[{3,4},"Double"],Bond[{5,6},"Double"],Bond[{6,7},"Single"],Bond[{2,7},"Single"]}}]&;
Grid[{{MoleculePlot[m1,IncludeHydrogens->False,ImageSize->144]," ⇌ ",MoleculePlot[m2,IncludeHydrogens->False,ImageSize->144]}},BaseStyle->{"Text",36,ShowStringCharacters->False}]
centerAndCloseCellGroup[]

The pattern reaction can be written conveniently with SMARTS:

\[KernelIcon];copeRearrangement=PatternReaction["[C:1]=[C:2][C:3][C:4][C:5]=[C:6]>>[C:3]=[C:2][C:1][C:6][C:5]=[C:4]"]

We can demonstrate this reaction with 1,5-heptadiene

The atoms have been numbered for convenience

\[KernelIcon];(heptadiene=Molecule["C=CCCC=CC"]//tagAtoms)//MoleculePlot

\[KernelIcon];ApplyReaction[copeRearrangement,{heptadiene}]//

The molecule bullvalene has three double bonds that are pairwise 1,5 to each other, and each pair can undergo a Cope rearrangement

\[KernelIcon];(bullvalene=Molecule["C2=CC1C=CC3C(C=C1)C23"])//

The arrangement of the 10 atoms is such that the rearrangement is degenerate, that is, the product is constitutionally the same as the reactant

This type of rearrangement is also known as tautomerization

coords={{-0.1966300000000003`,-0.35925`},{0.6939699999999993`,-0.35925`},{1.5801699999999999`,-0.013649999999999995`},{0.9597699999999998`,0.92995`},{-0.14783000000000035`,1.15595`},{-0.9438300000000002`,0.88075`},{-1.3500300000000003`,-0.25635`},{-0.5244300000000002`,-1.00165`},{0.6984699999999995`,-0.97505`},{-0.7696300000000003`,-0.001449999999999993`}};
m1=Molecule[ConstantArray["C",10],{Bond[{1,2},"Double"],Bond[{2,3},"Single"],Bond[{3,4},"Single"],Bond[{4,5},"Double"],Bond[{5,6},"Single"],Bond[{6,7},"Single"],Bond[{7,8},"Single"],Bond[{8,9},"Double"],Bond[{7,10},"Single"],Bond[{9,3},"Single"],Bond[{10,1},"Single"],Bond[{10,6},"Single"]},AtomDiagramCoordinates->coords];
(*MoleculePlot[m1,AtomLabels->Automatic]*)
m2=Molecule[ConstantArray["C",10],{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Single"],Bond[{4,5},"Double"],Bond[{5,6},"Single"],Bond[{6,7},"Single"],Bond[{7,8},"Double"],Bond[{8,9},"Single"],Bond[{2,9},"Single"],Bond[{9,3},"Single"],Bond[{10,1},"Double"],Bond[{10,6},"Single"]},AtomDiagramCoordinates->coords];
m3=Molecule[ConstantArray["C",10],{Bond[{1,2},"Single"],Bond[{2,3},"Single"],Bond[{3,4},"Single"],Bond[{4,5},"Single"],Bond[{5,6},"Double"],Bond[{6,7},"Single"],Bond[{7,8},"Single"],Bond[{8,9},"Double"],Bond[{7,10},"Single"],Bond[{9,3},"Single"],Bond[{10,1},"Double"],Bond[{4,2},"Single"]},AtomDiagramCoordinates->coords];
m4=Molecule[ConstantArray["C",10],{Bond[{1,2},"Double"],Bond[{2,3},"Single"],Bond[{3,4},"Single"],Bond[{4,5},"Single"],Bond[{5,6},"Double"],Bond[{4,9},"Single"],Bond[{7,8},"Double"],Bond[{8,9},"Single"],Bond[{7,10},"Single"],Bond[{9,3},"Single"],Bond[{10,1},"Single"],Bond[{10,6},"Single"]},AtomDiagramCoordinates->coords];
Column[{Row[{MoleculePlot[m1,IncludeHydrogens->False,ImageSize->144],MoleculePlot[m2,IncludeHydrogens->False,ImageSize->144]},Style["⇌",36]],Row[{MoleculePlot[m1,IncludeHydrogens->False,ImageSize->144],MoleculePlot[m3,IncludeHydrogens->False,ImageSize->144]},Style["⇌",36]],Row[{MoleculePlot[m1,IncludeHydrogens->False,ImageSize->144],MoleculePlot[m4,IncludeHydrogens->False,ImageSize->144]},Style["⇌",36]]}]
centerAndCloseCellGroup[]

The process exchanges all 10 carbon atoms; that is, there are 10!/3==1209600 isomers if all the atoms are labeled uniquely

\[KernelIcon];With[{mol=tagAtoms@bullvalene},
showReaction[mol,ApplyReaction[copeRearrangement,{tagAtoms@bullvalene},All]//Catenate,"Cope rearrangement",]
]

Dynamic NMR experiments by Oth et al. have confirmed the fluxional behavior of bullvalene, and that all the hydrogen atom positions are exchanged among themselves as are the carbon atom positions

At high temperature (>100°C) there is a single, average signal

As the temperature is lowered and the exchange process slows down, the signal broadens

At low temperature (~-60°C) distinct signals are observed for each chemically unique hydrogen or carbon position

(*Import["Oth74comparison.pdf","PageImages"][[4]]//ImageTake[#,900]&*)
Module[{cite},cite="Oth, J. F., Müllen, K., Gilles, J. M., & Schröder, G. (1974). Helvetica Chimica Acta, 57(5), 1415-1433.";
TabView[{
"(^1)H NMR" ->citation[Show[,ImageSize->12 72],cite],
"(^13)C NMR" ->citation[Show[,ImageSize->12 72],cite]},
ControlPlacement->Left
]
]
centerAndCloseCellGroup[]

We can study this process with token-event graphs by labeling one of the carbon atoms and observe it’s migration over all 4 symmetry nonequivalent positions

We’ll use a filled circle (●) in MoleculePlot for clarity

\[KernelIcon];(bullvaleneLabeledOne=Molecule["[13CH]2=CC1C=CC3C(C=C1)C23"])//moleculePlotLabeled

The chemically relevant token-event graph removes duplicate tokens (molecules) and events (reactions), and after 4 steps does not change

\[KernelIcon];g=TEG[{{b_}:>Sort@Catenate@ApplyReaction[copeRearrangement,{b},All]},{bullvaleneLabeledOne},4,"TokenDeduplication"->True,"EventDeduplication"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenLabeling"->True,"TokenRenderingFunction"->,ImageSize->8 72]

Labeling two positions demonstrates that all 12 symmetry nonequivalent pairs of positions can be labeled

\[KernelIcon];(bullvaleneLabeledTwo=Molecule["[13CH]2=[13CH]C1C=CC3C(C=C1)C23"])//moleculePlotLabeled

\[KernelIcon];g=TEG[{{b_}:>SortBy[#["IsomericSMILES"]&]@Catenate@ApplyReaction[copeRearrangement,{b},All]},{bullvaleneLabeledTwo},9,"TokenDeduplication"->True,"EventDeduplication"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenLabeling"->True,"TokenRenderingFunction"->,ImageSize->24 72,AspectRatio->1/2]

\[KernelIcon];Count[VertexList[g],_Molecule]

We can see that all possible pairs of positions have been labeled

\[KernelIcon];Cases[VertexList[g],_Molecule]//SortBy[#["IsomericSMILES"]&]//{Part[#,{3,2,1}],Part[#,{7,10,8,11,9}],Part[#,{5,12,6}],Part[#,{4}]}&//Map[moleculePlotLabeled,#,{2}]&//Grid
centerAndCloseCellGroup[]

Wagner-Meerwin Rearrangement

Carbocation Rearrangements

As alluded to in the introduction, the Wagner-Meerwein rearrangement involves carbocations, carbon atoms missing a bonded hydrogen atom and carrying a positive charge

Carbocations are also known as carbonium ions; there are four kinds:

MoleculePlot[#,ImageSize->144]&/@Molecule/@{"[CH3+]","[CH2+]CC(C)C","C[CH+]C(C)C","CC[C+](C)C"}//{#,{"methyl","primary or 1°","secondary or 2°","tertiary or 3°"}}&//Grid[#,Spacings->{2,2},BaseStyle->{"SmallText"}]&
centerAndCloseCellGroup[]

The stability of carbocations is: 3°>2°>1°>(SuperPlus[Me])

The order is due to the greater ability of alkyl groups than hydrogen atoms to donate elections toward the charged atom

Alkyl groups (methyl, ethyl, …) and hydrogen atoms vicinal to the cation tend to migrate, especially when a more stable cation is formed

This migration is also known as a 1,2-shift

1° carbocations are never formed by rearrangement

Here is an example using the iso-butyl cation

Module[{mp,m1,m2,m3,m4,eqa},
mp=MoleculePlot[#,IncludeHydrogens->"ExplicitOnly",ImageSize->144]&;
m1=Molecule[];
m2=Molecule[];
m3=Molecule[];
m4=Molecule[];
eqa=equilibriumArrow[{#,""},ImageSize->144,BaseStyle->{18}]&;
Row[{mp@m1,eqa["H migration"],mp@m2,eqa["Me migration"],mp@m3,eqa["H migration"],mp@m4}]
]
centerAndCloseCellGroup[]

Because Wagner-Meerwein rearrangements involve both carbon and hydrogen atoms and the propensity to form more stable carbocations, a single PatternReaction does not suffice and hydrogen atoms must be explicitly included

\[KernelIcon];(favoredAlkylShifts=PatternReaction/@<|"Alkyl shift 1->1"->,"Alkyl shift 1->2"->,"Alkyl shift 1->3"->,"Alkyl shift 2->2"->,"Alkyl shift 2->3"->,"Alkyl shift 3->2"->,"Alkyl shift 3->3"->|>)//

\[KernelIcon];(favoredHydrideShifts=PatternReaction/@<|"Hydride shift 1->1"->,"Hydride shift 1->2"->,"Hydride shift 1->3"->,"Hydride shift 2->2"->,"Hydride shift 2->3"->,"Hydride shift 3->2"->,"Hydride shift 3->3"->|>)//

Here is an example using the iso-pentane cation

The atoms have been numbered for convenience

\[KernelIcon];(isopentaneCation=Molecule["C[C+](C)CC"]//tagAtoms)//

\[KernelIcon];FoldList[First@ApplyReaction[#2,{#1}]&,isopentaneCation,{favoredHydrideShifts["Hydride shift 3->2"],favoredAlkylShifts["Alkyl shift 2->2"],favoredHydrideShifts["Hydride shift 2->3"]}]//

Norbornyl Cation Rearrangement

In the introduction we showed the rearrangement of borneol to camphene

Here, we will investigate the rearrangement of the norbornyl cation (bornane without the methyl groups)

The atoms have been numbered for convenience

\[KernelIcon];(norbornylCation=Molecule["C12[CH+]CC(C1)CC2"])//tagAtoms//

There are two possible 1,2-alkyl shifts: C:5 from C1: to C:2 and C:7 from C:1 to C:2

\[KernelIcon];With[{mol=tagAtoms@norbornylCation,key="Alkyl shift 2->2"},
showReaction[mol,ApplyReaction[favoredAlkylShifts[key],{mol},All]//Catenate,key,ImageSize->Small]
]
closeCellGroup[]

This is the first step of the token-event graph using  both the hydride- and alkyl-shift reactions, and without any deduplication

First, combine the PatternReactions into a single list of rules

\[KernelIcon];rules={cation_}:>SortBy[#["IsomericSMILES"]&]@Catenate@ApplyReaction[#,{cation},All]&/@Values@Join[favoredAlkylShifts,favoredHydrideShifts]//Echo[#,"|rules| =",Length]&;

We’ve colored the cation atom pink in the token rendering function to make it easier to follow

We used "EventLabeling"->"Colors" because there is more than one kind of reaction (rule)

\[KernelIcon];TEG[rules,{tagAtoms@norbornylCation},1,"EventDeduplication"->False,"EventLabeling"->"Colors","TokenDeduplication"->None,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenLabeling"->True,"TokenRenderingFunction"->,ImageSize->9 72]

From left to right:

the bond from C-7 to C-1 has shifted to C-2 and the charge shifted to C-1

the bond from C-5 to C-1 has shifted to C-2 and the charge shifted to C-1

a hydrogen atom on C-3 has shifted to C-2 and the charge shifted to C-3

a hydrogen atom on C-1 has shifted to C-2 and the charge shifted to C-1

In the absence of labeling, the first and third products are the same

The event coloring indicates which rule, and thus which reaction, was used

\[KernelIcon];(Keys@Join[favoredAlkylShifts,favoredHydrideShifts])[[{4,11,12}]]

Now, let’s construct a reaction network as a token-event graph

After 9 steps no new molecules are formed, and a total of 23 new molecules have been generated

For clarity we used "TokenDeduplication"->True and "EventDeduplication"->True

norbornylCationTEG=TEG[rules,{norbornylCation},9,"EventDeduplication"->True,"EventLabeling"->"Colors","TokenDeduplication"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenLabeling"->True,"TokenRenderingFunction"->,ImageSize->30 72]//Echo[#,"# molecules =",Count[VertexList@#,_Molecule]&]&

The reactions employed from one structure to another is easily interrogated from the graph, for example

\[KernelIcon];path=FindShortestPath[norbornylCationTEG,Molecule[],Molecule[]];

\[KernelIcon];showPath[path,Keys@Join[favoredAlkylShifts,favoredHydrideShifts],ImageSize->Tiny]

It’s somewhat easier to follow the evolution by suppressing the event vertices

TEG[rules,{norbornylCation},9,"EventDeduplication"->True,"TokenDeduplication"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"EventVertices"->False,"TokenLabeling"->True,"TokenRenderingFunction"->,ImageSize->10 72]//

The ring systems above the dashed line are found in naturally occurring monoterpenes

Notice that all but one of the steps are reversible, and it’s interesting to speculate that the role of the various terpene synthase enzymes may be to make the steps irreversible

Metabolism

Terpene Biosynthesis

In 1972 Devon and Scott published Handbook of Naturally Occurring Compounds

Volume II covers the terpenes, those molecules biosynthesized from isopentenyl pyrophosphate and dimethylallyl pyrophosphate via these precursors

Row[MapThread[Column[{MoleculePlot[#,ImageSize->#3],#2},Center]&,{(starToR[#,"PP"]&)@*Molecule/@{"CC(C)=CCO[*]","CC(C)=CCCC(C)=CCO[*]","CC(C)=CCCC(C)=CCCC(C)=CCO[*]","CC(C)=CCCC(C)=CCCC(C)=CCCC(C)=CCO[*]"},{"dimethylallyl PP","geranyl PP","farnesyl PP","geranylgeranyl PP"},72{1.75,2.25,3.25,3.75}}],DisplayForm[OverscriptBox["→",MoleculePlot[(starToR[#,"PP"]&)@Molecule["CC(=C)CCO[*]"],ImageSize->Tiny]]],BaseStyle->{"SmallText",18}]
centerAndCloseCellGroup[]

The figure shows the, at the time, best understood biosynthetic pathways of the various carbon skeletons

Module[{cite},cite="Devon, T. K. and Sott, A. I. (1972). Handbook of Naturally Occurring Compounds, Volume II, Terpenes. Academic Press.";
TabView[{
"Monoterpenes"->citation[Show[,ImageSize->12 72],cite],
"Sequiterpenes I"->citation[Show[,ImageSize->12 72],cite],
"Sequiterpenes II"->citation[Show[,ImageSize->12 72],cite],
"Diterpenes"->citation[Show[,ImageSize->12 72],cite]},
ControlPlacement->Left]
]
centerAndCloseCellGroup[]

We would like to use the methods covered here to fully elaborate the terpene biosynthetic landscape, and possibly even discover some unknown members

Strategy

One must be careful when writing the reactions

They cannot so promiscuous that unintended transformations occur

They need to avoid generating highly strained products

this is probably more easily handled as a post hoc filter on the output from ApplyReaction

They should closely mimic what the terpene synthase enzymes are known to do to avoid run-away isomer generation

The 1,2-shifts we studied above, along with with the reactions shown below, can be used to generate the transformations depicted above

{
{"cationic olefin alkylation",showReaction[Molecule["[CH2+]C=C(C)CCC=C(C)C"],Molecule["C1C=C(C)CCC1[C+](C)C"],ImageSize->Small]},
{"vinylogous cation rearrangement",showReaction[Molecule["[CH2+]C=C(C)[*]"]//starToR,Molecule["C=C[C+](C)[*]"]//starToR,ImageSize->Tiny]}
}//Grid[#,Alignment->{{Left,Center}},Spacings->3,BaseStyle->{"SmallText"}]&
centerAndCloseCellGroup[]

More PatternReactions

Cationic olefin alkylation

Note the use of recursive SMARTS to constrain the olefin to be at least two bonds away from a cation

\[KernelIcon];cationOlefinAlkylations=PatternReaction/@<|"C+ alkylation C=C"->"([C+:1].[C;!$(*[*+]):2]=[C;!$(*[*+]):3])>>([C+0:1][C:2][C+:3])"|>;

Demonstrate the reaction with geranyl cation

\[KernelIcon];With[{mol=Molecule["[CH2+]C=C(C)CCC=C(C)C"]},Map[showReaction[mol,#,ImageSize->Small]&,KeyValueMap[{#1, ApplyReaction[#2, mol, All]} &, cationOlefinAlkylations]//DeleteCases[{_,{}}],{3}]//MapAt[Column,#,{All,2}]&//Grid[#,Frame->All,FrameStyle->Gray,Spacings->{2,2}]&]
centerAndCloseCellGroup[]

Vinylogous cation rearrangement

\[KernelIcon];vinylogousCationRearrangements=PatternReaction/@<|
"VCR 1->1"->,"VCR 1->2"->,"VCR 1->3"->,
"VCR 2->2"->,"VCR 2->3"->,
"VCR 3->2"->,"VCR 3->3"->
|>;

Demonstrate the reaction with geranyl cation

\[KernelIcon];With[{mol=Molecule["[CH2+]C=C(C)CCC=C(C)C"]},Map[showReaction[mol,#,ImageSize->Small]&,KeyValueMap[{#1, ApplyReaction[#2, mol, All]} &, vinylogousCationRearrangements]//DeleteCases[{_,{}}],{3}]//MapAt[Column,#,{All,2}]&//Grid[#,Frame->All,FrameStyle->Gray,Spacings->{2,2}]&]
centerAndCloseCellGroup[]

Monoterpene Biosynthesis Token-Event Graph

\[KernelIcon];rules=({cation_}:>SortBy[#["IsomericSMILES"]&]@Catenate@ApplyReaction[#,{cation},All])&/@Values@Join[favoredAlkylShifts,favoredHydrideShifts,vinylogousCationRearrangements,cationOlefinAlkylations]//Echo[#,"|rules| =",Length]&;
ruleKeys=Keys@Join[favoredAlkylShifts,favoredHydrideShifts,vinylogousCationRearrangements,cationOlefinAlkylations];

\[KernelIcon];(geranylCation=Molecule["[CH2+]C=C(C)CCC=C(C)C"])//

Here are the first two steps, and already we can quite a bit of structural diversity

\[KernelIcon];TEG[rules,{geranylCation},2,"EventDeduplication"->True,"EventLabeling"->True,"TokenDeduplication"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenLabeling"->True,"TokenRenderingFunction"->,ImageSize->12 72]//Echo[#,"# molecules =",Count[VertexList@#,_Molecule]&]&

Using 6 steps returns a token-event graph with 1068 molecules, 76 of which (or their carbon skeleton) are present in the Devon and Scott monoterpenes shown above

I’ve taken the liberty to provide specific molecules in place of generic skeletons, and to replace some of the names with synonyms

With[{tegMols=Cases[VertexList[geranylCationTEGwithEvents],_Molecule]},Function[mol,MoleculePlot[mol,{Style[cationPosition[mol],Pink]},Background->Which[AnyTrue[tegMols,MoleculeMatchQ[mol,#]&],LightGreen,AnyTrue[tegMols,sameSkeletonQ[mol,#]&],LightYellow,True,Automatic],ImageSize->Tiny]/. {(*->Pink,*)White->Blue}]/@DevonAndScottMonoterpenes]//
centerAndCloseCellGroup[]

The molecules highlighted in green are present in the token-event graph, and those highlighted in yellow have a representative with the same carbon skeleton present in the token-event graph

The full token-event graph is too big to display usefully here, but we can trace the shortest paths from the root molecule geranyl cation to each of these 76 molecules

Module[{geranylCationTEGwithoutEvents,targets},
geranylCationTEGwithoutEvents=tokenGraph@geranylCationTEGwithEvents;
targets=Cases[VertexList[geranylCationTEGwithoutEvents],m_Molecule/;AnyTrue[DevonAndScottMonoterpenes,sameSkeletonQ[#,m]&]];
VertexDelete[geranylCationTEGwithoutEvents,Complement[VertexList[geranylCationTEGwithoutEvents],FindShortestPath[geranylCationTEGwithoutEvents, geranylCation, #] & /@ targets//Apply[Union]],GraphLayout->"LayeredDigraphEmbedding",ImageSize->24 72,AspectRatio->(1/2)]
]
closeCellGroup[]

How good is our assumption that our reactions mimic those used in nature?

The results are mixed:

iso-camphane skeleton

\[KernelIcon];isoCamphane=DevonAndScottMonoterpenes["iso-camphanyl cation"]//carbonSkeleton//#1["IsomericSMILES"]&//Echo[#1,"",]&

\[KernelIcon];targets=Select[VertexList[geranylCationTEGwithEvents],carbonSkeleton[#1]["IsomericSMILES"]==isoCamphane&]

\[KernelIcon];routes=(FindShortestPath[geranylCationTEGwithEvents,geranylCation,#1]&)/@targets//;

Most of the routes are extensions via hydride shifts of these two

\[KernelIcon];Map[showPath[#,ruleKeys,ImageSize->Tiny]&,routes[[{1,4}]]]//
centerAndCloseCellGroup[]

The first route is the accepted one (cf. Davis and Croteau)

The second route is doubtful

carane skeleton

\[KernelIcon];carane=DevonAndScottMonoterpenes["carene"]//carbonSkeleton//#1["IsomericSMILES"]&//Echo[#1,"",]&

\[KernelIcon];targets=Select[VertexList[geranylCationTEGwithEvents],carbonSkeleton[#1]["IsomericSMILES"]==carane&]

\[KernelIcon];routes=(FindShortestPath[geranylCationTEGwithEvents,geranylCation,#1]&)/@targets//;

The rest of the routes are extensions via hydride shifts of the first one

\[KernelIcon];Map[showPath[#,ruleKeys,ImageSize->Tiny]&,routes[[{1}]]]//
centerAndCloseCellGroup[]

Savage and Croteau have shown that pinyl cations are not intermediates:

citation[Show[,ImageSize->18 72],"Savage, T. J., & Croteau, R. (1993). Archives of biochemistry and biophysics, 305(2), 581-587."]
centerAndCloseCellGroup[]

Nucleophilic groups Y and Z are often invoked in terpene biosynthesis mechanisms, and are supposedly part of the enzymes in the synthase complex; and we have not included reactions utilizing them in this study

nepetalactol

DevonAndScottMonoterpenes["nepetalactol"]////
centerAndCloseCellGroup[]

This molecule’s skeleton was not discovered

Geu-Flores et al. have shown that the biosynthesis of iridoids proceeds via 8-oxogeranial and have proposed the following two mechanisms

citation[Show[,ImageSize->10 72],"Geu-Flores, F., Sherden, N. H., Courdavault, V., Burlat, V., Glenn, W. S., Wu, C., ... & O’Connor, S. E. (2012). Nature, 492(7427), 138-142."]
centerAndCloseCellGroup[]

Conclusions

Useful reactions have be written for the exploration of the monoterpene biosynthetic landscape

Reactions for the addition and expulsion of nucleophiles could be helpful

Additional reactions are needed for the classes of molecules not discovered

There is some prior art that we would be remiss if not mentioned:

Tian et al.

Studied monoterpenes

Used hydride- and alkyl-shifts, olefin alkylation, and proton transfer reactions

Geometry, stereochemistry, and energetics were considered

Exhaustive exploration

Tools: iGen (custom Python), Schrödinger API

da Silva et al.

Studied selected plant sesquiterpenes

Used selected olefin alkylations, deprotonation, [1,2]- and [1,3]-hydride shift, and hydration reactions

Used only graph rewriting

Limited exploration and ILP-based path discovery

Tools: MedØlDatschgerl, Cypher, Gremlin, CytoScape

Isegawa et al.

Studied sesquiterpenes

Reaction path investigation of humulyl cation using MO methods

Geometry, stereochemistry, and energetics of paramount importance

Tools: GRRM (Global Reaction Route Mapping), Gaussian09

Plans

Implement CC bond cleavage reactions that are known to be used, for example, for the chrysanthenamic acids

Study direct search (e.g., ResourceFunction["DepthFirstSearch"])and theorem proving methods to identify potential biosynthetic routes for specific input terpenes

Study the sesquiterpene and diterpene biosynthesis with the current method

Explore the use of stereospecific reactions

Acknowledgements

This presentation would not have been possible without the discussions, disagreements, support, and suggestions of the following people:

Wolfram Chemistry team

Jason Biggs

Jason Sonnenberg

Jan Mangaldan

Björn Zimmermann

Michael Trott

Wolfram Fundamental Physics Project team

Nik Murzin

Brad Klee

James Boyd

Stephen Wolfram

Thank you!

Initialization

Quit[]

SetDirectory@NotebookDirectory[];

closeCellGroup, centerAndCloseCellGroup

closeCellGroup[] := (SelectionMove[EvaluationCell[], Next,Cell]; 
FrontEndExecute[FrontEndToken["SelectionCloseUnselectedCells"]])

centerAndCloseCellGroup[] := (SelectionMove[EvaluationCell[], Next,Cell]; 
SetOptions[NotebookSelection[],{TextAlignment->Center,ShowCellLabel->False}];
FrontEndExecute[FrontEndToken["SelectionCloseUnselectedCells"]])

citation

citation[work_, str_String] := 
	Column[{
		work, 
		Style[str, "SmallText", FontSize -> 12, FontColor -> Gray]
		},
		Center, 
		Dividers -> {None,{False, Gray, False}}
		]

takeDropFold

takeDropFold[list_, nList:{___?NonNegative}] := 
	Fold[
		Module[{h, t}, {h, t} = TakeDrop[Last[#1], #2]; {Append[First[#1], h], t}] & , 
		{{}, list}, 
		nList
		] // Replace[#, {r:{_, {__}} :> Apply[Append, r], r_ :> First@r}]&

equilibriumArrow

ClearAll[equilibriumArrow]
Options[equilibriumArrow]={
"allowed"->True,
"preferred"->None, 
"style" -> Gray
};
equilibriumArrow[opts:OptionsPattern[{equilibriumArrow,Graphics}]]:=
equilibriumArrow[{"",""},opts]
equilibriumArrow[aboveLabel_String,opts:OptionsPattern[{equilibriumArrow,Graphics}]]:=
equilibriumArrow[{aboveLabel,""},opts]
equilibriumArrow[{aboveLabel_String,belowLabel_String},opts:OptionsPattern[{equilibriumArrow,Graphics}]]:=
equilibriumArrow[{aboveLabel,belowLabel},opts]
equilibriumArrow[{aboveLabel_String,belowLabel_String},opts:OptionsPattern[{equilibriumArrow,Graphics}]]:=
Catch@Module[{allowed,preferred},
allowed=OptionValue["allowed"];
If[!MatchQ[allowed,True|False],
Message[equilibriumArrow::opttf,"allowed",allowed];
Throw[$Failed]
];
preferred=OptionValue["preferred"];
If[!MatchQ[preferred,None|Left|Right],
Message[equilibriumArrow::optvg,"preferred",preferred,"None, Left, or Right"];
Throw[$Failed]
];
Graphics[
{
OptionValue["style"],
{Arrowheads[0.20],
If[preferred===Left,Arrow[{{-1.0,0.15},{1.0,0.15}}],Arrow[{{-1.5,0.15},{1.5,0.15}}]],
If[preferred===Right,Arrow[{{1.0,-0.15},{-1.0,-0.15}}],Arrow[{{1.5,-0.15},{-1.5,-0.15}}]],If[allowed,Nothing,Line[0.5{{{-1,1},{1,-1}},{{-1,-1},{1,1}}}]]},
{If[MatchQ[aboveLabel,Except[""]],Text[aboveLabel,{0,0.30},{0,-1}],Nothing],If[MatchQ[belowLabel,Except[""]],Text[belowLabel,{0,-0.30},{0,+1}],Nothing]}
},
FilterRules[{opts},Options@Graphics],
ImageSize->54,
PlotRangePadding->Scaled[0.1]
]
]

moleculeLayouts

moleculePlotLayouts[mol_Molecule,opts:OptionsPattern[MoleculePlot]]:=
Module[{methods={Automatic,"RDKit","CoordGen","OpenBabel","Avalon","PubChem"}},
Map[{ClickToCopy[MoleculePlot[Molecule[mol],AtomDiagramCoordinates->#,opts,ImageSize->Small],AtomDiagramCoordinates->#],#}&,{Automatic,"RDKit","CoordGen","OpenBabel","Avalon","PubChem"}]//Partition[#,UpTo[3]]&//Map[Transpose]//Apply[Join]//Grid[#,Dividers->{All,{{True,False},True}},ItemSize->Full]&
]

moleculePlotLayouts[Molecule["adamantane"]]

AtomDiagramCoordinates->"CoordGen"

starToR

starToR[mol_Molecule, label_:"R"] := Molecule[mol["AtomList"] /. 
    Atom["*"] -> Atom["*", "DisplayName" -> label], mol["BondList"]]

tagAtoms

Options@Symbol@First@Names["*`MoleculePlotMethod"]

"AtomMapNumber"

Options[tagAtoms]={IncludeHydrogens->False};
tagAtoms[mol_Molecule,opts:OptionsPattern[]]:=
Module[{labelHydrogensQ,patt,index,i=0},
labelHydrogensQ=TrueQ@OptionValue[IncludeHydrogens];
patt=Atom[If[labelHydrogensQ,_,Except["H"]],___];
index=Position[AtomList[mol,IncludeHydrogens->"ExplicitOnly"],patt];
Molecule[MapAt[Append[#,"AtomMapNumber"->++i]&,Replace[AtomList[mol],s_String:>Atom[s],{1}],index],BondList[mol],Options@mol]
]

tagAtoms[Molecule["alanine"]]//FullForm

tagAtoms[Molecule["alanine"]]//MoleculePlot[#,Method->{"ShowAtomMapNumbers"->True,"LabeledCarbons"->None}]&

"DisplayName"

Options[tagAtoms] = {IncludeHydrogens -> False}; 
tagAtoms[mol_Molecule, opts:OptionsPattern[]] := 
	Module[{labelHydrogensQ, patt, index, i = 0}, 
		labelHydrogensQ = TrueQ[OptionValue[IncludeHydrogens]]; 
		patt = Atom[If[labelHydrogensQ, _, Except["H"]], ___]; 
		index = Position[AtomList[mol, IncludeHydrogens -> "ExplicitOnly"], 
			patt]; 
		Molecule[
			MapAt[
				Append[#1, "DisplayName" -> StringTemplate["``:``"][First[#1], ++i]] & , 
				Replace[AtomList[mol], s_String :> Atom[s], {1}], 
				index
				], 
			BondList[mol], 
			Options[mol]
			]
		]

tagAtoms[Molecule["alanine"]]//FullForm

tagAtoms[Molecule["alanine"]]//MoleculePlot

tagAtoms[Molecule["alanine"],IncludeHydrogens->True]//MoleculePlot

cationPosition

cationPosition[mol_Molecule]:=FirstPosition[AtomList@mol,a_/;Lookup[Options@a,"FormalCharge",0]==1]

simpleLabel

simpleLabel[mol_Molecule] := 
	Molecule[
		mol["AtomList"] /. Atom[a_, opts1___, "MassNumber" -> m_, opts2___] :> 
			Atom[a, opts1, "MassNumber" -> m, "DisplayName" -> a, opts2], 
		mol["BondList"]
		]

moleculePlotLabeled

ClearAll[moleculePlotLabeled]
moleculePlotLabeled[mol_Molecule, opts:OptionsPattern[MoleculePlot]] := 
	MoleculePlot[
		Molecule[
			AtomList[mol] /. Atom[a_, "MassNumber" -> m_] :> 
				Atom[a, "DisplayName" -> "●"], 
			BondList[mol]
			], 
		opts, 
		AtomDiagramCoordinates -> "PubChem", 
		ImageSize -> Small
		]

showReaction

ClearAll[showReaction]
showReaction[m_Molecule, args___] := showReaction[{m}, args]
showReaction[r_, m_Molecule, args___] := showReaction[r, {m}, args]
showReaction[reactants:{__Molecule}, products:{__Molecule}, 
		agent_String:"         ", opts:OptionsPattern[MoleculePlot]] := 
	Row[Row[#, " + "]& /@ Map[MoleculePlot[#, opts]&, {reactants, products}, {2}], 
		DisplayForm[OverscriptBox["→", agent]], 
		BaseStyle -> {"SmallText", Gray}
		]

showPath

ClearAll[showPath]
showPath[path:{__}, ruleKeys:{__}, opts:OptionsPattern[MoleculePlot]] := 
	path /. {
		m_Molecule :> (MoleculePlot[m, {Style[cationPosition[m], Pink]}, 
			opts] /. {White -> Blue}), 
		{"Event", step_, ruleIndex_, __} :> 
			DisplayForm @ OverscriptBox["→", ruleKeys[[ruleIndex]]]
		} // Row

pubChemID, pubChemQ

ClearAll[pubChemID,pubChemQ]
pubChemID[mol_Molecule]:=pubChemID[mol]=mol["PubChemCompoundID"]
pubChemQ[mol_Molecule]:=pubChemQ[mol]=MatchQ[pubChemID[mol],{__}]
Scan[pubChemQ, ](* pre-cache for presentation *)

CASRegNo, CASRegQ

ClearAll[CASRegNo,CASRegQ]
CASRegNo[mol_Molecule]:=CASRegNo[mol]=mol["CASRegistryNumber"]
CASRegQ[mol_Molecule]:=CASRegQ[mol]=MatchQ[CASRegNo[mol],{__}]
Scan[CASRegQ, ](* pre-cache for presentation *)

carbonSkeleton, sameSkeletonQ

ClearAll[carbonSkeleton, sameSkeletonQ]

carbonSkeleton[molIn_String] := 
	With[{mol=Molecule[molIn]}, 
		carbonSkeleton[mol] /; MoleculeQ[mol]
		]

carbonSkeleton[molIn_Molecule] := carbonSkeleton[molIn] = 
	Module[{mol, atoms, bonds, p, b}, 
		atoms = AtomList[molIn]; 
		bonds = BondList[molIn]; 
		(* take only carbon atoms, discared all properties *)
		p = First /@ Position[atoms, Atom["C", ___], {1}]; 
		atoms = atoms[[p]] /. Atom[a_, ___] :> Atom[a]; 
		(* take only bonds between carbon atoms, set type to single and renumber *)
		b = Thread[p -> Range@Length[p]]; 
		bonds = Select[bonds, AllTrue[First@#, MemberQ[p, #]&]&] /. Bond[idx_, ___] :> Bond[idx /. b, "Single"]; 
		mol = Molecule[atoms, bonds]; 
		(* take the largest fragment *)
		ConnectedMoleculeComponents[mol] // TakeLargestBy[AtomCount, 1] // First
		]

sameSkeletonQ[mols__Molecule] := sameSkeletonQ[mols] = 
	carbonSkeleton /@ {mols} // AllTrue[Partition[#, 2, 1], Apply[MoleculeMatchQ]]&
sameSkeletonQ[] := False

sameSkeletonQ[Molecule["cyclohexane"],Molecule["cyclohexene"],Molecule["benzene"]]

sameSkeletonQ[Molecule[],Molecule[]]

TEG

TokenEventGraph is under active development, and is available from three sources:

nikTEG=ResourceFunction[CloudObject["https://www.wolframcloud.com/obj/nikm/\
DeployedResources/Function/TokenEventGraph"]];(* described (HyperlinkTemplate[here, https://www.wolframcloud.com/obj/nikm/DeployedResources/Function/TokenEventGraph/, https://www.wolframcloud.com/obj/nikm/DeployedResources/Function/TokenEventGraph/, HyperlinkActionRecycled, HyperlinkActive, , ]) *)
wppTEG=ResourceFunction["https://www.wolframcloud.com/obj/wolframphysics/\
DeployedResources/Function/TokenEventGraph"]; (* described (HyperlinkTemplate[here, https://www.wolframcloud.com/obj/wolframphysics/DeployedResources//Function/TokenEventGraph, https://www.wolframcloud.com/obj/wolframphysics/DeployedResources//Function/TokenEventGraph, HyperlinkActionRecycled, HyperlinkActive, , ]) *)
wfrTEG=ResourceFunction["TokenEventGraph"]; (* described (HyperlinkTemplate[here, https://resources.wolframcloud.com/FunctionRepository/resources/TokenEventGraph/?i=TokenEventGraph&searchapi=https%3A%2F%2Fresources.wolframcloud.com%2FFunctionRepository%2Fsearch, https://resources.wolframcloud.com/FunctionRepository/resources/TokenEventGraph/?i=TokenEventGraph&searchapi=https%3A%2F%2Fresources.wolframcloud.com%2FFunctionRepository%2Fsearch, HyperlinkActionRecycled, HyperlinkActive, , ]) *)
TEG = Switch["WPP",
	"Nik", nikTEG, (* takes list of states, which are lists of tokens *)
	"WPP", wppTEG, (* takes list of tokens *)
	_(* WFR *), wfrTEG (* takes list of tokens *)
]

Nik’s options

FilterRules[Options@TEG,Except[Options@Graph]]//Column

WPP’s options

FilterRules[Options@TEG,Except[Options@Graph]]//Column

wppTEG[{{diene_,dienophile_}:>Catenate@ApplyReaction[dielsAlder,{diene,dienophile},All]},{cyclohexadiene,isoprene},1,"TokenDeduplication"->True,"TokenLabeling"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenRenderingFunction"->(MoleculePlot[#,ImageSize->72]&),ImageSize->12 72]

nikTEG[{{diene_,dienophile_}:>Catenate@ApplyReaction[dielsAlder,{diene,dienophile},All]},{{cyclohexadiene,isoprene}},1,"TokenDeduplication"->True,"TokenLabeling"->True,"TokenEquivalenceFunction"->MoleculeMatchQ,"TokenRenderingFunction"->(MoleculePlot[#,ImageSize->72]&),ImageSize->12 72]

takeStepsTokenEventGraph

takeStepsTokenEventGraph[teg_Graph, maxStep_Integer?NonNegative, opts:OptionsPattern[Graph]] := 
	Module[{v, e, graphLayout, performanceGoal, vertexLabels, vertexShapeFunction, vertexSize, vertexStyle}, 
		{v, e} = If[maxStep == 0, 
			FirstCase[EdgeList[teg], DirectedEdge[token_, {"Event", 0, __}, ___] :> {{token}, {}}, {{}, {}}], 
		  (* else *)
			e = Cases[EdgeList[teg], DirectedEdge[_, {"Event", s_ /; s < maxStep, __}, ___] | 
				DirectedEdge[{"Event", s_ /; s < maxStep, __}, _, ___]]; 
			{Union @@ (Take[#, 2]&) /@ List @@@ e, e}
			]; 
		
		{graphLayout, performanceGoal, vertexLabels, vertexShapeFunction, vertexSize, vertexStyle} = 
			{GraphLayout, PerformanceGoal, VertexLabels, VertexShapeFunction, VertexSize, VertexStyle} /. Options[teg]; 
		vertexLabels = Cases[vertexLabels, HoldPattern[(k_ -> _) /; MemberQ[v, k]]]; 
		vertexShapeFunction = Cases[vertexShapeFunction, HoldPattern[(k_ -> _) /; MemberQ[v, k]]]; 
		vertexSize = Cases[vertexSize, HoldPattern[(k_ -> _) /; MemberQ[v, k]]]; 
		vertexStyle = Cases[vertexStyle, HoldPattern[(k_ -> _) /; MemberQ[v, k]]]; 
		
		Graph[v, e, opts, 
			GraphLayout -> graphLayout, PerformanceGoal -> performanceGoal, 
			VertexLabels -> vertexLabels, VertexShapeFunction -> vertexShapeFunction, 
			VertexSize -> vertexSize, VertexStyle -> vertexStyle
			]
		]

tokenGraph

remove event vertices from the graph

useful when events have only one input

no checking is done for this condition

ClearAll[tokenGraph]
Options[tokenGraph]={
"untagEdges"->False
};
tokenGraph[g_Graph,opts:OptionsPattern[{tokenGraph,Graph}]]:=
Module[{v,e,gOpts=Association@Options@g},
v=Cases[VertexList@g,_Molecule];
e=Join@@(Cases[EdgeList@g,DirectedEdge[#[[2]],token_,tag_]:>DirectedEdge[#[[1]],token,tag]]&)/@Cases[EdgeList@g,DirectedEdge[_Molecule,__]];
If[TrueQ@OptionValue["untagEdges"],
e=e/.DirectedEdge[a_,b_,_]:>DirectedEdge[a,b]//DeleteDuplicates
];
gOpts[VertexLabels]=DeleteCases[gOpts[VertexLabels],{"Event",__}->_];
gOpts[VertexStyle]=DeleteCases[gOpts[VertexLabels],{"Event",__}->_];
Graph[v,e,Sequence@@FilterRules[{opts},Options[Graph]],Sequence@@Normal@gOpts]
]

extractRouteDAR

ClearAll[extractRouteDAR]
extractRouteDAR[l:{___Molecule}, args___] := 
	(extractRouteDAR[#1, args] & ) /@ l
extractRouteDAR[mol_Molecule, opts:OptionsPattern[TEG]] := 
	Module[{teg, paths, edges}, 
		teg = TEG[{{diene_, dienophile_} :> Catenate[ApplyReaction[dielsAlder, 
			{diene, dienophile}, All]]}, {}, 2, 
			"TokenDecompositionFunction" -> 
				(DeleteDuplicates[Join[#1, {cyclohexadiene, isoprene}], MoleculeMatchQ] & ), "EventDeduplication" -> True, 
			"TokenDeduplication" -> True, 
			"TokenEquivalenceFunction" -> MoleculeMatchQ, 
			"TokenLabeling" -> True, 
			"TokenRenderingFunction" -> 
				(MoleculePlot[#1, AtomDiagramCoordinates -> "CoordGen", 
					ImageSize -> 72] & ), 
			opts]; 
		paths = Catenate[(FindPath[teg, #1, mol, Infinity, Infinity] & ) /@ 
			{cyclohexadiene, isoprene}]; 
		edges = (Partition[#1, 2, 1] & ) /@ paths // Catenate // DeleteDuplicates // 
			MapApply[Cases[EdgeList[teg], DirectedEdge[##1, _]] & ] // Catenate; 
		Subgraph[teg, edges]
		]

norbornylCationTEG

norbornylCationTEG = ;

VertexCount[norbornylCationTEG]
EdgeCount[norbornylCationTEG]

DevonAndScottMonoterpenes

(DevonAndScottMonoterpenes=Molecule/@<|
(* row 1 *)
"cyclolavandulal"->"C1CC(C)(C)CC(C)=C1(C=O)",
"santolina triene"->"C=C(C)CC(=C(C)C)C=C",
"chrysanthemic acid"->"CC(C)=CC1C(C1(C)C)C(=O)O",
"artemesia ketone"->"CC(C)=CC(=O)C(C)(C)C=C",
(* row 2 *)
"beta-pyronene"->"C1C=CC(C)=C(C)C1(C)(C)",
"thujyl cation"->"C1C[C+](C)C2CC12C(C)C",
"lavandulol"->"CC(C)=CCC(CO)C(C)=C",
"carene"->"C1CC(C)=CC2C(C)(C)C12",
"eucarvone"->"C1=CC=C(C)C(=O)CC1(C)(C)",
(* row 3 *)
"eucarvone "->"C1=CC=C(C)C(=O)CC1(C)(C)",
"geranyl cation"->"[CH2+]C=C(C)CCC=C(C)C",
"α-terpinyl cation"->"C1CC(C)=CCC1[C+](C)C",
"tropone"->"C1=CC(=O)C=CC(C(C)C)=C1",
"sylvecarvone"->"C1CC(=O)C(C)=CC1C(C)=C",
(* row 4 *)
"nepetalactol"->"C12CCC(C)C1C(O)OC=C2C",
"bornyl cation"->"C12CCC(C)(C1(C)(C))[CH+]C2",
"pinyl cation"->"C12CC[C+](C)C(C1(C)(C))C2",
"fenchyl cation"->"C12CCC(C)(C1)[CH+]C2(C)(C)",
(* row 5 *)
"seco-iridoid"->"CCC(C(C)C)C(C)CC",
"iso-camphanyl cation"->"C12CCC(C1)[C+](C)C2(C)(C)",
"teresantolol"->"C12CC3C(C1)C3(C)C2(C)(CO)",
"carquejol"->"C1C=CC(=C)C(C(C)=C)C1O",
(* row 6 *)
"1-acyl-4-propenyl-1-cyclopentene"->"C1C(C(=O)C)CC(C(C)=C)=C1",
"filfolone"->"C12CC=C(C)C1C(=O)C2(C)(C)"
|>)//Grid[MapThread[Insert[#1,"",#2]&,{takeDropFold[KeyValueMap[Column[{#2,#1},Center]&,MoleculePlot[#,ImageSize->Tiny]&/@#1],{4,5,5,4,4,2}],{{{1}},{},{},{{1}},{{1}},{{1},{2},{2}}}}],BaseStyle->"SmallText"]&

geranylCationTEGwithEvents

Module[{stream},
stream=OpenRead["geranylCationTEG_withEvents_2--2022-10-07-20-51-35.wxf"];
geranylCationTEGwithEvents=ReadList[stream]//Last;
Close[stream];
Echo@Count[VertexList[geranylCationTEGwithEvents],_Molecule];
]

$Line

$Line=0;