# Chemistrees:
data driven identification of reaction pathways via machine learning.

Sander Roet  [<EMAIL_ADDRESS>Christopher David Daub  [  Enrico Riccardi
[

###### Abstract

We propose a supervised machine learning algorithm, decision trees, to analyze
molecular dynamics output. The approach aims to identify the predominant
geometric features which correlate with trajectories that transition between
two arbitrarily defined states. The data-based algorithm aims to identify such
features in an approach which is unbiased by human ”chemical intuition”.

We demonstrate the method by analysing proton exchange reactions in formic
acid (FA) solvated in small water clusters. The simulations were performed
with _ab initio_ molecular dynamics combined with a method for generating rare
events, specifically path sampling. Our machine learning analysis identified
mechanistic descriptions of the proton transfer reaction for the different
water clusters.

###### keywords:

Chemometrics, Decision Tree, Classifiers, Molecular Dynamics, Path Sampling

NTNU]Department of Chemistry, Norwegian University of Science and Technology,
Trondheim, Norway UofH]Department of Chemistry, University of Helsinki, P.O.
Box 55, FI-00014, Helsinki, Finland UiO]Department of Informatics, UiO,
Gaustadalléen 23B, 0373 Oslo, Norway

##  1 Introduction

In regions far from urban areas, formic acid (FA) has been recognised as one
of the main factors which reduces the pH of rainwater, causing acid rain. 1 It
has relatively high atmospheric concentrations 2, 3, and contributes to the
formation of sulphuric acid in the atmosphere 4, 5. Enhanced description and
prediction of proton exchange reactions involving solvated FA can improve
current atmospheric models.

_Ab initio_ molecular dynamics simulations have recently been used to examine
FA deprotonation in aqueous solution 1, 6, demonstrating that simulations can
successfully model proton exchange between water and FA. While these studies
led to valuable new insights, the limitations of the adopted methods (_e.g._
usage of a bias potential, continuous collective variables) could be overcome
thanks to relatively novel methodologies, such as Replica Exchange Transition
Interface Sampling (RETIS) 7, 8. Respecting the natural dynamics of the system
allows the study of transitions even with a significant diffusive contribution
9, 10 (_i.e._ a small reaction barrier), and enables the direct investigation
of reaction mechanisms.

RETIS is a rare event method developed to investigate transitions. Its main
advantages are: (a) it does not alter the natural dynamics of the system, (b)
it does not require a particularly accurate order parameter, (c) its results
are in principle identical to what would be obtained by an infinitely long
unbiased molecular dynamics simulations. With RETIS, the transition region is
explored by continuously generating new paths which start from a stable state
and end up either back in such a state – an _unreactive_ path, or reach a
different state – a _reactive_ path. The approach has been successfully
employed to study transitions that would, otherwise, require prohibitively
long simulation times. The results generated have been used to describe the
dynamics of chemical processes (_e.g._ reaction rates) while considering the
entropic contribution in the analysis11, 12, 6, 9, 10. Since large quantities
of data are often generated by the sampling procedure, approaches to
pragmatically decode reaction mechanisms are greatly beneficial.

Our aim is to establish an heuristic approach to describe transitions,
regardless of whether they involve crossing an entropic barrier. Data driven,
physically consistent and measurable system descriptors might be generated,
and their correlations with the system dynamics asserted. This is a
classification problem which a Machine Learning (ML) algorithm can be trained
to solve. The algorithm might then predict if a certain molecular structure
(frame) is part of a reactive or a non-reactive trajectory. Connecting the
descriptors to measurable quantities provides a data driven (”unbiased”)
description of a transition that might support, and eventually surpass, human-
biased ”chemical intuition”.

Data driven algorithms for enhanced sampling or the analysis of chemical
simulations have significant recent contributions. 13, 14, 15, 16, 17, 18.
Most of these approaches are based on neural networks, which lack physically
consistent interpretability. By contrast, this can be a desirable feature of
decision trees 15, 16. The usage of neural networks in rare events,
furthermore, requires a pre-selection of trial collective variables 14, 17,
18, which could lead to a hypothesis bias. Decision tree 19 (DT) classifiers
have an unique solution and are not sensitive to highly correlated variables.
The results can be readily interpreted if the source variables are also
interpretable. This approach was previously adopted to select optimal
collective variables with decision trees, with reasonable success. 12

Our work based on decision trees provides both an interpretable and
hypothesis-bias-free method and proposes the usage of interpretable DTs via an
appropriate system representation invariant to system translation, rotation,
and changes in atomic indices. Our aim is to gain insights into reaction
mechanisms with a systematic and unbiased representation of the system.

The approach has been developed with sufficient versatility to be applied to
different types of molecular simulations, from conventional molecular dynamics
(MD) to rare event methods. It shall be noted that conventional MD would
require _a priori_ classification of the data (_i.e._ dividing the source
trajectory into reactive and unreactive segments). The sampling strategy of
rare event methods, instead, generates a data structure which inherently
classifies the trajectories. Regardless of the adopted molecular simulation
approach, limiting the correlation between samples is a primary task for a
quantitative data driven method to identify reaction paths and the probability
of their occurrence.

In the present study, we focus on small clusters of formic acid solvated by
water, HCOOH + (H2O)n, $n=4$ and $6$. The system is relatively small and well
understood, and hence provides an ideal test case for training a ML method. At
the same time, studying the system using our methodology provides new
quantitative and qualitative insights into the acid-water proton transfer
reaction which is also of interest in atmospheric chemistry. 20, 21, 22, 23,
24, 25, 26.

##  2 Computational Models and Methods

Our main focus is the machine learning methodology here introduced, and
therefore only a brief introduction to the simulation methodology will be
provided. Please consider our previous work6, 5 for further details.

###  2.1 System description

For studying proton transport, molecular simulations able to consider bond
formation and bond breaking are required. Born-Oppenheimer molecular dynamics
(BOMD) has been shown to be a suitable approximation in previous studies of
atmospheric reactions 4, 5 and of aqueous formic acid 1, 6. The density
functional theory BLYP implemented in the Quickstep module of CP2K 27 has been
combined with a double-zeta basis set supplemented by the use of Grimme’s D2
dispersion correction. 28

Figure 1: Minimum energy configurations for systems with formic acid (FA)
associated with 4 and 6 water molecules. $s(r)$ in all cases shown here is
equivalent to the $r_{OH}$ of the initially protonated formic acid (FA)
molecule. These configurations correspond to the initial states used to start
PyRETIS simulations.

A set of systems with an increasing number of water molecules around FA were
studied. Initial configurations were obtained from minimum energy
configurations, of which snapshots are reported in Figure 1. As more water
molecules were added, the probability of generating reactive trajectories
became larger. However, at least four added water molecules were required to
allow generation of trajectories with a significant charge separation between
the deprotonated FA and the solvated proton. Adding two more water molecules
resulted in a significantly higher proton transfer rate. The two systems
composed by the formic acid surrounded by four and six water molecules are
here analysed and discussed.

###  2.2 Definition of the Collective Variable

Path sampling simulations require the definition of a collective variable,
$s(r)$, where $r$ represents the set of positions and velocities of all atoms
in the system. The collective variable must describe the transition from a
reactant state A to a product state B. Path sampling is not limited to
continuous variables, which allows the consideration of relatively more
complex functions to describe proton transport.

The collective variable here adopted, $s(r)$, is inspired by the study of
water ionization 12, with modifications introduced to consider acid
deprotonation. As a first step, it locates the smallest distance between any
formic acid oxygen and any reactive hydrogen in the system excluding the
methyl hydrogen in formic acid. This distance is denoted as $r_{O_{FA}H,min}$.

For $r_{O_{FA}H,min}<1.4$ Å the formic acid can be considered protonated, so
$s(r)=r_{O_{FA}H,min}$. For $r_{O_{FA}H,min}>1.4$ Å, charge separation between
the solvated proton and the formic acid becomes significant. To quantify it,
and thus compute $s(r)$, all the distances between reactive hydrogens and
oxygens are first considered. Hydrogens are then assigned to the closest
oxygen, either water or formic acid. Any water oxygen found to be associated
with 3 hydrogens is indexed. Next, all distances between formic acid oxygens
and hydrogens associated with triply coordinated water oxygens are sorted.
$s(r)$ is the minimum value of these distances.

Conceptually, the reasoning behind this complex definition of $s(r)$ is to
allow the description of formation of more extended systems such as Eigen or
Zundel cations. A discontinuous jump of $s(r)$ from $\sim 1.8$ Å up to $\sim
3$ Å is associated with a change in the identity of the triply coordinated
water oxygen and the formation of structures resembling Zundel cations
H5O${}_{2}^{+}$. The formation of the Zundel cation with $s(r)>2.9$ Å is here
labeled as the product state B.

###  2.3 RETIS

The PyRETIS 29, 30 library has been used to perform RETIS31 simulation,
coupled with the ab-initio molecular dynamics external engine, CP2K. The
initial paths describing the transition from protonated to deprotonated FA
along $s(r)$ were generated by using the _kick_ method available in the
software, starting from the initial configurations shown in Figure 1. The
_kick_ approach uses a mixture of stochastic and deterministic dynamics to
generate a set of initial paths. The results were obtained by discarding paths
correlated with the initial generated ones. Finally, the remaining paths from
three and four independent simulations for the four and six water cases ,
respectively, were merged together.

The RETIS simulations varied somewhat in the details of the number of
interfaces used for the transition interface sampling, and the location of the
first interface (_i.e._ the interface between the $[0^{-}]$ and the $[0^{+}]$
ensemble). In the 4 water simulations, the first interface was placed at
$s(r)=1.05$ Å, and the last interface at $s(r)=3.0$ Å. Similarly, the 6 water
simulations had the first interface at $s(r)=1.07$ Å, and the last interface
at $s(r)=3.0$ Å.

###  2.4 Selection window

In a trajectory, each frame can be considered as an instance. Depending on the
simulations set up, a large number of frames would generate a long list of
instances with a very high correlation. Furthermore, different trajectories
can be highly correlated one another, depending on the sampling algorithm.
Since generating a sufficient number of uncorrelated trajectories often
requires excessive computational requirements, an approach to provide a
sufficient sampling with a limited correlation is here proposed. We have thus
selected frames contained in a rather restricted region in path space, the
selection windows. By randomly picking a certain number of frames, for each
trajectory, within a selection window, the correlation between instances is
minimized. Depending on the location and size of the selection window,
different outcomes can be expected. By placing it, as in the current study, in
proximity of the initial state, the system configurations that are correlated
with the transitions can be identified at an early stage of the transition.
The selection window location and dimension variables, along with the number
of frames per trajectory to consider, constitute the three hyper-parameters of
our approach. In the present work, the machine learning algorithm has been
thus fed with one frame per trajectory with order parameters $1.1<s(r)<1.25$
Å. The range is sufficiently narrow to consider only a few frames for each
trajectory, each with similar order parameter.

The selected frames, each of them belonging to a different trajectory, have
comparable values of the order parameter. The machine learning algorithm
should, hereby, determine the most relevant feature(s) for the transition to
happen, without bias on the main descriptor of the transition itself. This
ensures that the detected features are not correlated to the classification of
the trajectory.

###  2.5 Training of the Decision trees

Decision trees report the most important features to differentiate between
reactive and unreactive paths, without imposing any prior hypothesis. In the
present work, we considered atomic distances and velocities as possible
features. Both entries have to be treated such to be translation and
rotational invariant. As we did not achieve any significant improvement in the
decision tree by considering atomic velocities, we here report the work based
on the analysis of atomic distance only.

In most of the simulation packages, Molecular simulations produce an ordered
data array for the positions and velocities for each atom in a time frame.
While the convention facilitates post processing and visualization procedures,
it includes a bias in the data representation. Small deviations in the
observation angle, or on the axis definition (e.g. exchanging $x$ with $y$
coordinates) leads to significantly different data sets, while corresponding
to nearly identical systems. Data thus has to be pre-processed to become
invariant with respect to translations and rotations. Furthermore, the machine
learning problem also has to be atom-index invariant, and the sorting method
also must be reversible to allow back-mapping of the machine learning
prediction to the relevant atom (or atom pairs).

The translation and rotation independent requisites are met by a simple atom-
atom distance matrix. For the FA acid in a water cluster the carbon atom is
the trivial identifiable reference atom. The atom-index invariant requisite
is, on the other hand, a more elaborate representation. The matrix is created
by first choosing an anchor point (C0) and grouping the row per element and
sorting the distance with the anchor atom for each element. The resulting
matrix reports the distance from a selected anchor atom (rows) to its next
neighbour (columns). A schematic of the algorithm to generate both the simple
distance matrix and the atom invariant distance matrix is provided in the SI.
In Figure 2, the two matrix representations are provided, as an example, for a
formic acid molecule.

| C0 | H0 | O0 | O1 | H1
---|---|---|---|---|---
C0 | 0.0000 | 0.1109 | 0.1250 | 0.1321 | 0.2010
H0 | 0.1109 | 0.0000 | 0.2047 | 0.2013 | 0.2954
O0 | 0.1250 | 0.2047 | 0.0000 | 0.2311 | 0.2524
O1 | 0.1321 | 0.2013 | 0.2311 | 0.0000 | 0.1073
H1 | 0.2010 | 0.2954 | 0.2524 | 0.1073 | 0.0000

| C0 | H0 | H1 | O0 | O1
---|---|---|---|---|---
C0 | 0.0000 | 0.1109 | 0.2010 | 0.1250 | 0.1321
H0 | 0.1109 | 0.0000 | 0.2954 | 0.2013 | 0.2047
H1 | 0.2010 | 0.0000 | 0.2954 | 0.1073 | 0.2524
O0 | 0.1250 | 0.2047 | 0.2524 | 0.0000 | 0.2311
O1 | 0.1321 | 0.1073 | 0.2013 | 0.0000 | 0.2311

Figure 2: Distance matrix (top) for a structure of formic acid (right). The
index invariant distance matrix (bottom) corresponding to the first distance
matrix.

The resulting internal coordinate representation allows an independent
analysis of each entry and, thus, a suitable data structure for the machine
learning task.

Another common translational and rotational invariant representation, the
Z-matrix 32, also provides an appealing internal representation of molecular
structures as it scales better with the number of atoms compared to the
distance matrix. Yet, the value of its variables (_i.e._ distances, angles and
dihedrals) are dependent on each other and on the atom sequence. In contrast,
in our distance matrix representation, the considered variables are all
independent. Also, distances are unique, with a lower (0) bound and an upper
bound (system size). These three characteristics allow for a proper split of
sample space to be made by the decision trees. Furthermore, our representation
can be index invariant.

The machine learning problem we are posing is: what are the main features that
a simulation frame has to have in order to be part of a trajectory that
connects an initial state to a final state (reactive)? With which probability?

The information gain (entropy) decision tree is a viable method for a problem
with highly correlated features. 33

Given a set of trajectories, a classification between reactive and unreactive
paths is first needed. A numerical descriptor, conventionally defined as order
parameter, can quantify the progress of a given transition. If its value, for
a given system, is within certain arbitrarily defined ranges, the system can
be considered to be located in the initial or final state. A reactive path is
defined as a path starting from an initial state and ending at a product
state. A non reactive path, instead, ends at the initial state. If the input
generated by molecular simulation is composed by a single long trajectory,
sub-segments will have to be fed to the machine learning task. If a segment
starting at one state and ending in another state will be considered reactive,
while a segment starting and ending at the same state without having entered
another state before, is considered as unreactive. When using the input
generated by path sampling, paths contained in a single ensemble should be
considered (please consider Refs. 7, 34 for the definition of an ensemble and
further details of the path sampling methods). Regardless of the approach to
generate input paths, or segments, it is worth noting here that the population
of reactive and unreactive paths has to be sufficient to feed a machine
learning algorithm. The necessary population size depends largely on the
system size and type. In the attempt to provide an estimate of the reliability
of the predictions, an error-estimate procedure has been also developed and
described in a forthcoming section.

Computationally, a DecisionTree Classifier from scikit-learn 35 has been fed
with the index invariant matrix, using the ’entropy’ splitting criterion and a
max depth of 3.

###  2.6 Visualizing decision paths

A symmetric distance matrix can be back-mapped to $xyz$ coordinates (up to a
translation and rotation), as described by Young and Householder 36 (and
further detailed in the SI). The index-invariant distance matrix can be
unsorted into the symmetric matrix, up to an atom index difference. The
approach permits to add dummy atoms according to the splits given by the
decision tree, allowing a direct visualization of the analysis output (e.g.
via VMD 37). For a convenient visualization, only the dummy atoms
corresponding to the nodes along each decision path in the tree are selected.
A main decision path is chosen such that the leaf node would have the highest
number of pertinent reactive paths, weighted by the percentage of pertinent
reactive paths: $n_{r}\cdot\frac{n_{r}}{n_{r}+n_{u}}$, where $n_{r}$ is the
number of reactive paths in that node and $n_{u}$ the number of unreactive
paths.

###  2.7 Random forest decision error-estimate

With highly correlated data, significantly different trees can be originated
depending upon the first split from minor variations of the input. This is a
constitutive limitation of the approach. The relative importance of the first
split can be asserted by using Random Forests with a unit depth. The feature
importance of such Random Forests is equal to the importance of the first
split only. IT shall be noted that the feature that has the highest importance
in the Random forest plot does not necessarily represent the main split for
all possible decision trees.

A sequence of forests of decision trees have been generated in respect with
the time sequence on which sampling output has been generated. Source data has
been split in 10 sub-blocks and randomized within. Random forest has been then
computed for each of these subsets, generating a sort of time-dependent
profile for the main splits, which allows to compute a variance $\sigma$ on
each of the main features. The average value, for each feature, can be instead
computed by considering the whole data set. By assuming a Gaussian
distribution for each feature, and by using the previously obtained variance
and mean value for each feature, the relative probability of a feature
importance can be estimated. By comparing the probability distribution for
each feature, the most relevant can be identified even in the presence of
highly correlated data.

Computationally, a RandomForest Classifier from scikit-learn 35 has been fed
with the index invariant matrix, using the ’entropy’ splitting criterion and a
max depth of 1.

##  3 Results and discussion

The rate of proton transfer from the formic acid to the water molecules has
been computed via RETIS simulation and ab-initio molecular dynamics
simulations. Figure 3 reports the rate of reaction for two systems, where 4
and 6 water molecules surrounded the FA acid. State A (protonated state) is
defined as configurations with $s(r)<1.05$ Å(4 waters) or with $s(r)<1.07$ Å(6
waters). State B includes configurations with $s(r)>3.0$ Åfor both systems.

Figure 3: Effective rate constant $k_{AB}$ computed using RETIS for FA
clustered with 4 or 6 molecules to reach the deprotonated state. Results are
an average of several RETIS simulations from different initial conditions,
weighted by the number of RETIS cycles in each simulation.

The rate of proton transfer for the four water molecule case is
$~{}2.1010^{-}141/fs$, while $~{}1.0110^{-}71/fs$ for six water molecules
around the formic acid ($~{}10^{7}$ times difference).

We here investigate the mechanism of reactions via decision trees to identify
the feature(s) that better correlate, for each case, with the proton transfer.
The analysis might provide qualitative and quantitative description of the
different system features responsible for the significant difference in the
reported rates.

For the two systems, we performed different independent simulations using
different initialization procedures available within the PyRETIS package 30 to
reduce the correlation in the exploration of path space and to speed up
sampling. The results obtained from each independent simulation have been then
combined. Even if the analysis of individual simulations further demonstrated
the efficiency of our method to reveal partial difference within realizations
(e.g. different, less likely, reaction channels), we limit ourselves to
discussion of the findings from the total combined data in the present
manuscript.

Generally, all trajectory segments, or trajectories for path sampling, can be
considered in the present analysis approach. When using the latter simulation
approach, a re-weight algorithm shall be adopted to consider all the generated
paths. Due to the statistical weight of the different ensemble, and for
simplicity, we opted to consider only the trajectories included in the the
outermost ensemble in path space (for the definition of an ensemble please
consider the RETIS formalism 31).

Each considered trajectory, or trajectory segment, has to be labeled according
to a reactive/unreactive definition. The frames from the selection windows
have been transformed into the distance matrix representation (transnational
and rotational invariant) and fed to the machine learning algorithm to
generate the forthcoming analysis. We want to stress here that these results
are generated without any arbitrary hypothesis or bias. The only hyper-
parameters are the number of frames per trajectory, the location and size of
the selection window.

Decision trees, visualisations and random forest of decision trees to assert
the first split error have been generated for the two systems. We here report
the results obtained by the index-invariant distance matrix, which is the most
general approach, even if more computationally demanding. It is worth noting
that the index-variant distance matrix can be still advantageous for its
simplicity and symmetry in certain applications. For example, in the presence
of atoms that do now swap order during a transition. In the SI, the results
for the index-variant distance matrix will also be presented.

The atom labelling system we use identifies each atom with a character and a
digit. The character corresponds to the atom type, while the digit corresponds
to the position of the sorted distance list per element with respect to a
reference atom, with the indexing starting at 0.

The digit of the first entry in the atom-atom distance label refers to the
sorted distance list with respect to the reference atom (C). The digit of the
second entry refers to the sorted distance list with respect to the first atom
of the atom-atom pair. As an example, O2-H5 corresponds to the distance of the
third closest oxygen (O2) from the C atom to the hydrogen atom which is 6th
closest to the oxygen in the pair (O2). H0-O0 is the distance from the H
closest to the C (H0) to the oxygen closest to H0.

###  3.1 Formic Acid with 4 water molecules

The decision tree generated for the four water molecules cluster around the FA
is reported in figure 4. To simplify the visualization of the main splits that
lead to the highest reactive trajectories of the decision tree, in the figure
4, a xyz representation has been included. The atoms in blue, yellow, and
green are involved in the first, second, and third split, respectively.

The deprotonation reaction of FA appears to primarily require the distance
between O5 and H9 being smaller than 5.25 Å. This split implies that the
distance between the furthest oxygen and the furthest hydrogen from such
oxygen should be within a given threshold. As H9 is the hydrogen of FA, that
also implies that a certain orientation of the molecule, in respect to the
water cluster, is also required. In these conditions, the probability for the
path to be reactive is 38%.

The next split, along the branch which leads to the highest reactive stance,
is the distance between O1 and H8 being smaller than 4.25 Å. The distance
between one of the FA oxygens and one of the furthest hydrogen should be
sufficiently small. That implies that the oxygen of the FA should be located
around the center of the cluster and that a sort of specific structure of the
water cluster is required to satisfy this condition. With this condition, the
probability for a path to be reactive reaches 63%

Still along the branch towards the highest reactive stance, the distance
between O5 and 03 being bigger than 3.52 Å represent the last split here
considered. This corresponds to the relative position of two water molecules
being two hydrogen bonds away. We interpret this requisite as the suitable
distance to establish hydrogen bonding between the atoms. When the three
requirements are met, the probability of for a path of being reactive is of
71%

By comparing the number of reactive paths versus the number of unreactive
paths in the final splits of the decision tree, it can be concluded that the
indicated reactive path is clearly predominant. A similar conclusion can be
reached by observing the first splits reported in figure 5. The figure, which
reports the results obtained from a random forest of decision trees of depth
1, indicates the relative probability for the first split to be the most
important feature with a 39% probability. The following distances reported by
the random forest have a constantly decaying relevance. The first 5 main
splits reported by figure 5 are correlated, and indicate that the water
cluster has to be sufficiently compact and the FA has to be oriented such that
its oxygen are in contact with the surrounding water molecules.

Figure 4: Decision tree for the system with 4 water molecules around the
formic acid molecule based on the index-invariant distance matrix. Three
splits divide the input, each square reports 1) the question to split the
data, 2) the number of samples going into the node, 3) the number of
[unreactive, reactive] samples going into the node and 4) the majority class
of the node. At each split, the True branch is on the left, False on the
right. The color indicates the ratio between unreactive (brown) and reactive
(blue) samples included in a node. Wider arrows have been used to link the
decision three split with the atoms involved. In the top left corner, a 3D
representation of the system is provided. The atoms highlighted in blue,
yellow, and green, correspond to the atoms involved in the first, second, and
third split of the decision three, respectively. In red, white and back are
the oxygen, hydrogen and carbon atoms not indicated by the decision tree.
Figure 5: The importance approximations of the first split question from a
random forest with depth of one for the four water case system. The bars
represent the feature importance of a Random Forest, with the error bar
calculated with a block-error average based on the generated trajectories. The
probability that a split is the true most important split is shown above the
bars for the 3 most probable first splits.

###  3.2 Formic Acid with 6 water molecules

For the six water molecules water clusters around FA, in figure 6, we report
the generated decision tree. Consistently with the four water molecule case, a
visualization of the main splits of the decision tree that lead to the highest
reactive trajectories is also included in figure 6. The atoms in blue, yellow,
and green are involved in the first, second, and third slit, respectively.

The deprotonation reaction of FA in the six water molecule cluster requires a
distance between O6 and H8 being smaller than 3.55 Å. This split involves two
water molecules in the proximity of FA, that needs to be within a certain
distance. By inspecting the frame reported in figure 6, the requirement seems
to indicate a water structure in the surrounding of the water molecule in
proximity with the FA oxygen. In such conditions the probability for the path
to be reactive is 32%.

The next split, along the branch that leads to the highest reactive stance, is
the distance between H9 and H11 being smaller than 4.27 Å. The distance
between these two atoms can also be interpreted as a combination of molecular
orientation of the water molecules in the surrounding of the FA and water
cluster size. The probability of a reactive path reaches 53% when such
conditions occur.

Still along the branch towards the highest reactive stance, the distance
between O2 and 02 being bigger than 2.57 Å represent the last split here
considered. This indicates that the closest water oxygen to FA (O2) should be
close enough to a second oxygens to promote the formation of a hydrogen bond
network. When the three requirements are met, the probability for a path of
being reactive is of 72%

By comparing the number of reactive paths versus the number of unreactive
paths in the final splits of the decision tree, it can be concluded that the
indicated reactive path is clearly favourable, but that other significant
paths can also exist. This conclusion is also supported by the random forest
of decision trees with a single split. Before proposing an interpretation, it
is worth to remind here that the random forest reports unconditional entries,
while the decision tree splits depending from the first split. Figure 7
indicates the first split to be the most important feature is 34%, while the
second has a comparable relevance: O2-O2 (28%). This confirms that, while a
predominant pathway for the reaction has been sampled, different main pathways
can co-exist.

Figure 6: Decision tree for the system with 6 water molecules around the
formic acid molecule. Three splits divide the input, each square reports 1)
the question to split the data, 2) the number of samples going into the node,
3) the number of [unreactive, reactive] samples going into the node and 4) the
majority class of the node. At each split, the True branch is on the left,
False on the right. The color indicates the ratio between unreactive (brown)
and reactive (blue) samples included in a node. Wider arrows have been used to
link the decision three split with the atoms involved. In the top left corner,
a 3D representation of the system is provided. The atoms highlighted in blue,
yellow, and green, correspond to the atoms involved in the first, second, and
third split of the decision three, respectively. In red, white and back are
the oxygen, hydrogen and carbon atoms not indicated by the decision tree.
Figure 7: The importance approximations of the first split question from a
random forest with depth of one for the 6 water case system. The bars
represent the feature importance of a Random Forest, with the error bar
calculated with a block-error average based on the generated trajectories. The
probability that a split is the true most important split is shown above the
bars for the 3 most likely first splits.

As reported in figure 3, the dimension and the water cluster had a significant
effect on the rate of the proton transfer reaction. From the comparison of the
previously discussed figures 4 and 6, it can be noted that the distance
between an oxygen and a hydrogen in the furthermost regions in respect to the
FA is the predominant characteristic for a trajectory to be reactive. Both
clusters have to be sufficiently compact in order to promote the reaction. In
the four water molecule case, the orientation of the FA in respect of the
water cluster is the most important feature, while for the six water molecule
case the water structure around the FA appears to be the predominant feature.

A second main difference between the four and six water cases is the possible
pathways for the reaction to occur. The smaller system has a predominant
reactive path, while for the six water molecule cluster, multiple paths appear
to co-exist, contributing to the final reaction rate. Physically, if the
system is sufficiently large, different configurations can lead to the proton
transfer reaction.

##  4 Conclusions

A data driven method to systematically compute reaction pathways has been
presented. The conventional xyz data representation employed in molecular
simulations is converted into an index invariant distance matrix
representation, which is also translation and rotation invariant. Thereafter,
an approach to limit the correlation between source data (molecular dynamic
trajectories) has been also proposed, in consistency with rare event
simulation framework. The data has then fed to a supervised classifier method,
the decision tree. To simplify the interpretation of the classifier, a back
mapping procedure from the index invariant matrix has been adopted to
emphasize the atoms involved with each split identified by the decision tree.
Random forest of decision trees, in combination with block averaging, provided
an error range of the first split of the decision tree.

We thus presented a data-driven approach to gain insight into a chemical
reaction. The method has been designed such to be readily applicable to other
simulation strategies and transitions. The strength of the present approach is
to generate complex collective variables, and estimate the probability of
their occurrence in a transition path. The descriptors to elucidate transition
mechanisms might be directly implemented in a prediction method. 13

The method adopted an index invariant distance matrix providing a data driven
insight into the reaction pathways. The data driven identification aims to
identify interpretable pathways in a system composed of indistinguishable
molecules. Applications to more inhomogeneous systems would be
straightforward, especially if only a portion of the system atom is of
interest. The latter case would combine a human intuition with a data driven
approach, which would, possibly, provide a better insight into the reaction
if, and only if, the introduced bias is correct. Our method can be further
expanded by considering a higher number of descriptors alongside the distance
matrix. Velocities, angles between molecules, coarse grain procedure or a mix
of user defined functions, 38 can be fed to the decision tree and follow up
analysis.

To demonstrate the capabilities of the developed method, a mechanistic
description of the proton transfer reaction has been provided. The reaction
has been simulated via rare event simulation (replica exchange transition
interface sampling 7) and its rate quantified for two water clusters, one
composed of four and one of six water molecules surrounding a FA molecule.

The reaction rate we measure is strongly influenced by the number of water
molecules present. Mechanistically, the four and six water proton transfer
reaction requires the water cluster to be sufficiently compact. The four water
molecules system requires a certain orientation of the FA molecules and of the
water molecules in its proximity. For the six water molecule case, a certain
orientation of the outer water molecules appears to be more significant in
describing the reaction path. Furthermore, the four water cluster system
indicated a predominant pathway for the reaction to occur, while in the six
water molecules cluster, several pathways have been identified contributing to
the reaction.

Part of the work has been performed under the Project HPC-EUROPA3
(INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action
under the H2020 Programme; in particular, the authors gratefully acknowledge
the support of the Department of Chemistry at the University of Helsinki and
the computer resources and technical support provided by CSC. Christopher D.
Daub acknowledges funding by the Academy of Finland (Grant Number 294752) and
by the Jane and Aatos Erkko Foundation.

## References

  * Murdachaew et al. 2016 Murdachaew, G.; Nathanson, G. M.; Gerber, R. B.; Halonen, L. Deprotonation of formic acid in collisions with a liquid water surface studied by molecular dynamics and metadynamics simulations. _Phys. Chem. Chem. Phys._ 2016, _18_ , 29756–29770
  * Millet et al. 2015 Millet, D. B. et al. A large and ubiquitous source of atmospheric formic acid. _Atmos. Chem. Phys._ 2015, _15_ , 6283–6304
  * Chaliyakunnel et al. 2016 Chaliyakunnel, S.; Millet, D. B.; Wells, K. C.; Cady-Pereira, K. E.; Shephard, M. W. A Large Underestimate of Formic Acid from Tropical Fires: Constraints from Space-Borne Measurements. _Environ. Sci. Technol._ 2016, _50_ , 5631–5640
  * Kangas et al. 2020 Kangas, P.; Hänninen, V.; Halonen, L. An Ab Initio Molecular Dynamics Study of the Hydrolysis Reaction of Sulfur Trioxide Catalyzed by a Formic Acid or Water Molecule. _J. Phys. Chem. A_ 2020,
  * Daub et al. 2020 Daub, C. D.; Riccardi, E.; Hänninen, V.; Halonen, L. Path sampling for atmospheric reactions: formic acid catalysed conversion of SO3 $+$ H2O to H2SO4. _PeerJ Physical Chemistry_ 2020, _2_ , e7
  * Daub and Halonen 2019 Daub, C. D.; Halonen, L. Ab Initio Molecular Dynamics Simulations of the Influence of Lithium Bromide Salt on the Deprotonation of Formic Acid in Aqueous Solution. _J. Phys. Chem. B_ 2019, _123_ , 6823–6829
  * van Erp 2007 van Erp, T. S. Reaction Rate Calculation by Parallel Path Swapping. _Phys. Rev. Lett._ 2007, _98_ , 268301
  * Riccardi et al. 2017 Riccardi, E.; Dahlen, O.; van Erp, T. S. Fast Decorrelating Monte Carlo Moves for Efficient Path Sampling. _J. Phys. Chem. Lett._ 2017, _8_ , 4456–4460
  * Riccardi et al. 2019 Riccardi, E.; Van Mastbergen, E. C.; Navarre, W. W.; Vreede, J. Predicting the mechanism and rate of H-NS binding to AT-rich DNA. _PLoS computational biology_ 2019, _15_ , e1006845
  * Riccardi et al. 2020 Riccardi, E.; Krämer, A.; van Erp, T. S.; Ghysels, A. Permeation Rates of Oxygen through a Lipid Bilayer Using Replica Exchange Transition Interface Sampling. _The Journal of Physical Chemistry B_ 2020,
  * Moqadam et al. 2017 Moqadam, M.; Riccardi, E.; Trinh, T. T.; Lervik, A.; van Erp, T. S. Rare event simulations reveal subtle key steps in aqueous silicate condensation. _Physical Chemistry Chemical Physics_ 2017, _19_ , 13361–13371
  * Moqadam et al. 2018 Moqadam, M.; Lervik, A.; Riccardi, E.; Venkatraman, V.; Alsberg, B. K.; van Erp, T. S. Local initiation conditions for water autoionization. _Proc. Nat. Acad. Sci._ 2018, _115_ , E4569–E4576
  * van Erp et al. 2016 van Erp, T. S.; Moqadam, M.; Riccardi, E.; Lervik, A. Analyzing complex reaction mechanisms using path sampling. _Journal of Chemical Theory and Computation_ 2016, _12_ , 5398–5410
  * Hooft et al. 2021 Hooft, F.; de Alba Ortíz, A. P.; Ensing, B. Discovering Collective Variables of Molecular Transitions via Genetic Algorithms and Neural Networks. _Journal of Chemical Theory and Computation_ 2021,
  * Chen and Ferguson 2018 Chen, W.; Ferguson, A. L. Molecular enhanced sampling with autoencoders: On-the-fly collective variable discovery and accelerated free energy landscape exploration. _Journal of Computational Chemistry_ 2018, _39_ , 2079–2102
  * Schöberl et al. 2019 Schöberl, M.; Zabaras, N.; Koutsourelakis, P.-S. Predictive collective variable discovery with deep Bayesian models. _The Journal of Chemical Physics_ 2019, _150_ , 024109
  * Ribeiro et al. 2018 Ribeiro, J. M. L.; Bravo, P.; Wang, Y.; Tiwary, P. Reweighted autoencoded variational Bayes for enhanced sampling (RAVE). _The Journal of Chemical Physics_ 2018, _149_ , 072301
  * Jung et al. 2019 Jung, H.; Covino, R.; Hummer, G. Artificial intelligence assists discovery of reaction coordinates and mechanisms from molecular dynamics simulations. _arXiv preprint arXiv:1901.04595_ 2019,
  * Swain and Hauska 1977 Swain, P. H.; Hauska, H. The decision tree classifier: Design and potential. _IEEE Transactions on Geoscience Electronics_ 1977, _15_ , 142–147
  * Leopold 2011 Leopold, K. R. Hydrated Acid Clusters. _Annu. Rev. Phys. Chem._ 2011, _62_ , 327–349
  * Forbert et al. 2011 Forbert, H.; Masia, M.; Kaczmarek-Kedziera, A.; Nair, N. N.; Marx, D. Aggregation-Induced Chemical Reactions: Acid Dissociation in Growing Water Clusters. _J. Am. Chem. Soc._ 2011, _133_ , 4062–4072
  * Chung and Kim 2017 Chung, Y. K.; Kim, S. K. Dissociation of sulfur oxoacids by two water molecules studied using ab initio and density functional theory calculations. _Int. J. Quantum Chem._ 2017, _117_ , e25419
  * Lengyel et al. 2017 Lengyel, J.; Pysanenko, A.; Fárník, M. Electron-induced chemistry in microhydrated sulfuric acid clusters. _Atmos. Chem. Phys._ 2017, _17_ , 14171–14180
  * Gutberlet et al. 2009 Gutberlet, A.; Schwaab, G.; Birer, Ö.; Masia, M.; Kaczmarek, A.; Forbert, H.; Havenith, M.; Marx, D. Aggregation-Induced Dissociation of HCl(H2O)4 Below 1 K: The Smallest Droplet of Acid. _Science_ 2009, _324_ , 1545
  * Maity 2013 Maity, D. K. How Much Water Is Needed To Ionize Formic Acid? _J. Phys. Chem. A_ 2013, _117_ , 8660–8670
  * Elena et al. 2013 Elena, A. M.; Meloni, S.; Ciccotti, G. Equilibrium and Rate Constants, and Reaction Mechanism of the HF Dissociation in the HF(H2O)7 Cluster by ab Initio Rare Event Simulations. _J. Phys. Chem. A_ 2013, _117_ , 13039–13050
  * VandeVondele et al. 2005 VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. QUICKSTEP: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. _Comput. Phys. Commun._ 2005, _167_ , 103–128
  * Grimme et al. 2010 Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. _The Journal of chemical physics_ 2010, _132_ , 154104
  * Lervik et al. 2017 Lervik, A.; Riccardi, E.; van Erp, T. S. PyRETIS: A well-done, medium-sized python library for rare events. _Journal of Computational Chemistry_ 2017, _38_ , 2439–2451
  * Riccardi et al. 2020 Riccardi, E.; Lervik, A.; Roet, S.; Aarøen, O.; van Erp, T. S. PyRETIS 2: an improbability drive for rare events. _Journal of computational chemistry_ 2020, _41_ , 370–377
  * van Erp et al. 2003 van Erp, T. S.; Moroni, D.; Bolhuis, P. G. A novel path sampling method for the calculation of rate constants. _J. Chem. Phys._ 2003, _118_ , 7762–7774
  * Baker and Hehre 1991 Baker, J.; Hehre, W. J. Geometry optimization in cartesian coordinates: The end of the Z-matrix? _Journal of computational chemistry_ 1991, _12_ , 606–610
  * Breiman et al. 2017 Breiman, L.; Friedman, J. H.; Olshen, R. A.; Stone, C. J. _Classification And Regression Trees_ ; Routledge, 2017
  * Cabriolu et al. 2017 Cabriolu, R.; Refsnes, K. M. S.; Bolhuis, P. G.; van Erp, T. S. Foundations and latest advances in replica exchange transition interface sampling. _J. Chem. Phys._ 2017, _147_ , 152722
  * Pedregosa et al. 2011 Pedregosa, F. et al. Scikit-learn: Machine Learning in Python. _Journal of Machine Learning Research_ 2011, _12_ , 2825–2830
  * Young and Householder 1938 Young, G.; Householder, A. S. Discussion of a set of points in terms of their mutual distances. _Psychometrika_ 1938, _3_ , 19–22
  * Humphrey et al. 1996 Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. _J. Molec. Graphics_ 1996, _14_ , 33–38
  * Aarøen et al. 2021 Aarøen, O.; Kiær, H.; Riccardi, E. PyVisA: Visualization and Analysis of path sampling trajectories. _Journal of Computational Chemistry_ 2021, _42_ , 435–446