\makeatletter \@ifundefined{HCode} {\documentclass[screen,CRCHIM,Unicode,Thematic]{cedram} \newenvironment{noXML}{}{} \def\tsup#1{$^{{#1}}$} \def\tsub#1{$_{{#1}}$} \def\thead{\noalign{\relax}\hline} \def\tbody{\noalign{\relax}} \def\endthead{\noalign{\relax}\hline} \def\ndash{\text{--}} \RequirePackage{etoolbox} \usepackage{upgreek} \def\jobid{crchim20240075} %\graphicspath{{/tmp/\jobid_figs/web/}} \graphicspath{{./figures/}} \newcounter{runlevel} \let\MakeYrStrItalic\relax \csdef{Seqnsplit}{\\} \def\refinput#1{} \def\back#1{} \def\0{\phantom{0}} \newenvironment{inftab}{}{} \usepackage{multirow} \def\nrow#1{\@tempcnta #1\relax% \advance\@tempcnta by 1\relax% \xdef\lenrow{\the\@tempcnta}} \def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}} \DOI{10.5802/crchim.317} \datereceived{2024-01-31} \daterevised{2024-05-14} \dateaccepted{2024-05-15} \ItHasTeXPublished } {\documentclass[crchim]{article} \def\bibliographystyle#1{} \def\bibliography#1{} \let\refinput\input \def\CDRdoi{10.5802/crchim.317} \def\0{\phantom{0}} \let\newline\break \makeatletter } \makeatother \usepackage{upgreek} \dateposted{2024-11-27} \begin{document} \newcommand{\molgraph}{2D-MolGraph\xspace} \newcommand{\zala}{Z-Ala$_6$--COOH\xspace} \newcommand{\chondro}{Chondroitin disulfate CS2S4S\xspace} \newcommand{\xchondro}{chondroitin disulfate CS2S4S\xspace} \begin{noXML} %\makeatletter %\def\TITREspecial{\relax} %\def\cdr@specialtitle@english{French Network on Solvation (GDR 2035 SolvATE)} %\def\cdr@specialtitle@french{R\'eseau th\'ematique sur la solvatation (GDR 2035 SolvATE)} %\makeatother \title{Topological graphs: a review of some of our achievements and perspectives in physical chemistry and homogeneous catalysis} \alttitle{Graphes topologiques : une revue de nos r\'{e}alisations et des perspectives en chimie physique et en catalyse homog\`{e}ne} \author{\firstname{Sana} \lastname{Bougueroua}\CDRorcid{0000-0002-0642-5275}\IsCorresp} \address{Universit\'e Paris-Saclay, Univ Evry, CY Cergy Paris Universit\'e, CNRS, LAMBE UMR8587, 91025 Evry-Courcouronnes, France} \email[S. Bougueroua]{sana.bougueroua@univ-evry.fr} \author{\firstname{Yl\`ene} \lastname{Aboulfath}} \address{Universit\'e Paris-Saclay, Univ Versailles Saint Quentin, DAVID, 78035 Versailles, France} \author{\firstname{Alvaro} \lastname{Cimas}\CDRorcid{0000-0001-8095-2279}} \addressSameAs{1}{Universit\'e Paris-Saclay, Univ Evry, CY Cergy Paris Universit\'e, CNRS, LAMBE UMR8587, 91025 Evry-Courcouronnes, France} \author{\firstname{Ali} \lastname{Hashemi}} \address{Inorganic Systems Engineering, Department of Chemical Engineering, Faculty of Applied Sciences, Delft University of Technology, 2629 HZ Delft, The Netherlands} \author{\firstname{Evgeny~A.} \lastname{~Pidko}\CDRorcid{0000-0001-9242-9901}} \addressSameAs{3}{Inorganic Systems Engineering, Department of Chemical Engineering, Faculty of Applied Sciences, Delft University of Technology, 2629 HZ Delft, The Netherlands} \author{\firstname{Dominique} \lastname{Barth}} \addressSameAs{2}{Universit\'e Paris-Saclay, Univ Versailles Saint Quentin, DAVID, 78035 Versailles, France} \author{\firstname{Marie-Pierre} \lastname{Gaigeot}\CDRorcid{0000-0002-3409-5824}\IsCorresp} \addressSameAs{1}{Universit\'e Paris-Saclay, Univ Evry, CY Cergy Paris Universit\'e, CNRS, LAMBE UMR8587, 91025 Evry-Courcouronnes, France} \address{Institut Universitaire de France (IUF), 75005 Paris, France} \email[M.-P. Gaigeot]{mgaigeot@univ-evry.fr} \shortrunauthors \keywords{\kwd{Algorithmic graph theory} \kwd{Conformational search} \kwd{Molecular dynamics} \kwd{Identification of conformers} \kwd{Pathways} \kwd{Reaction network}} \altkeywords{\kwd{Th\'{e}orie algorithmique des graphes} \kwd{Recherche conformationnelle} \kwd{Dynamiques mol\'{e}culaires} \kwd{Identification de conform\`{e}res/conformations} \kwd{Chemins} \kwd{R\'{e}seau de r\'{e}activit\'{e}}} \begin{abstract} This paper reviews some of our developments in algorithmic graph theory, with some applications in physical chemistry and catalysis. Two levels of granularity in the topological graphs have been developed: atomistic 2D-MolGraphs and coarse-grained polygraphs of H-bonded cycles. These graphs have been implemented with the key algorithms of isomorphism and polymorphism, in order to analyze molecular dynamics simulations of complex molecular systems. These topological graphs are transferable without modification from ``simple'' gas molecules, to liquids, to more complex inhomogeneous interfaces between solid and liquid for instance. We show hereby that the use of algorithmic graph theory provides a direct and fast approach to identify the actual conformations sampled over time in a trajectory. Graphs of transitions can also be extracted, showing at first glance all the interconversions over time between these conformations. H-bond networks in condensed matter molecular systems such as aqueous interfaces are shown to be easily captured through the topological graphs. We also show how the 2D-MolGraphs can easily be included in automated high-throughput \mbox{in silico} reactivity workflows, and how essential they are in some of the decisive steps to be taken in these workflows. The coarse-grained polygraphs of H-bonded cycles are shown to be essential topological graphs to analyze the dynamics of flexible molecules such as a hexapeptide in\unskip\break gas phase. \end{abstract} \begin{altabstract} Cet article passe en revue certains de nos d\'{e}veloppements dans le domaine de la th\'{e}orie algorithmique des graphes, avec quelques applications en physico-chimie et en catalyse. Deux niveaux de granularit\'{e} dans les graphes topologiques ont \'{e}t\'{e} d\'{e}velopp\'{e}s : les 2D-MolGraphs atomistiques et les polygraphes \`{a} gros grain des cycles \`{a} liaisons hydrog\`{e}nes. Ces graphes ont \'{e}t\'{e} mis en \oe uvre avec les algorithmes cl\'{e}s de l'isomorphisme et du polymorphisme, pour l'analyse et l'interpr\'{e}tation de trajectoires de dynamique mol\'{e}culaire de syst\`{e}mes mol\'{e}culaires complexes. Ces graphes topologiques peuvent \^{e}tre transf\'{e}r\'{e}s sans modification des mol\'{e}cules pr\'{e}sentes en phase gazeuse aux liquides, en passant par des interfaces inhomog\`{e}nes plus complexes entre un solide et un liquide, par exemple. Nous montrons ici que l'utilisation de la th\'{e}orie algorithmique des graphes fournit une approche directe et rapide pour identifier les conformations r\'{e}elles qui sont \'{e}chantillonn\'{e}es au fil du temps dans une trajectoire. Des graphes de transitions peuvent \'{e}galement \^{e}tre extraits, montrant d'un seul coup d'\oe il toutes les interconversions au fil du temps entre ces conformations. Les r\'{e}seaux de liaisons hydrog\`{e}nes dans les syst\`{e}mes mol\'{e}culaires de la mati\`{e}re condens\'{e}e, tels que les interfaces aqueuses, sont facilement captur\'{e}s par les graphes topologiques. Nous montrons \'{e}galement comment les 2D-MolGraphs peuvent \^{e}tre facilement inclus dans les flux de travail automatis\'{e}s de r\'{e}activit\'{e} in silico \`{a} haut d\'{e}bit, et \`{a} quel point ils sont essentiels pour certaines des \'{e}tapes d\'{e}cisives \`{a} franchir dans ces flux de travail. Les polygraphes \`{a} gros grain des cycles \`{a} liaisons hydrog\`{e}nes se r\'{e}v\`{e}lent \^{e}tre des graphes topologiques essentiels pour l'analyse de la dynamique des mol\'{e}cules flexibles telles qu'un hexapeptide en phase gazeuse. \end{altabstract} \maketitle \vspace*{1.8pc} \twocolumngrid \end{noXML} \let\rmmu\upmu \section{Introduction}\label{sec1} The field of Theoretical and Computational Chemistry applies the laws of physics and chemis\-try coupled with computer programs to calculate the structures and chemical and physical properties of molecules in different states of matter, such as thermodynamic properties, spectroscopic signals, chemical reaction pathways, phase diagrams, etc. The arsenal of theoretical tools in computational chemistry has evolved in recent years, nowadays also including theoretical methods from Operation Research (OR), which uses algorithms to build solutions on well-formulated problems, or Artificial Intelligence (AI), which uses various Machine Learning methods, based in particular on neural networks, not only to predict new physical and chemical states, events and properties, but also to develop, for example, force fields or DFT functionals for simulations. This new era has triggered a revival in the field of theoretical and computational chemistry for research teams to develop new theoretical methods that include OR and/or AI to go beyond the simple use of ``classical numerical methods''. Over the past decade, our group has approached this new era through the prism of algorithmic graph theory, in the context of OR and AI, based on the representation of matter in topological graphs~\cite{bougueroua2018graph,bougueroua2023algorithmic,MolPhys2023,JCTC_ylene_2024,hashemi2022renegate}. A graph encodes topological properties of matter (in the same way for molecules, assemblies of molecules, liquids, solid materials, also interfaced with liquids) by means of vertices and edges that reflect the specific interactions (in pairs) between vertices. At the molecular level of representation, the vertices are usually associated to atoms while the edges report on interactions between atoms, e.g., chemical bonds and intermolecular interactions. Most graphs are defined in two dimensions (2D-graphs), any information related to, e.g., distances, angles, is usually not encoded into topological 2D-graphs unless vertices/nodes are specifically labeled with such information. Graphs can however be also three-dimensional, with an indication of coordinates in space that would thus encode intra- and inter-molecular interactions between atoms. Easier to obtain or to predict than 3D-graphs, 2D-graphs already carry information on the structure, the functional properties, and even the 3D shape of the materials they model. Examples include the\unskip\break classification of similar molecules according to their topology~\cite{NW19,GD13}, the prediction of patterns in biological molecules~\cite{BV15}, the prediction of the 3D structure of small molecules~\cite{LD13}, etc. Molecular graphs are also commonly used in supervised machine learning algorithms, the framework of graph-based models for molecules is indeed naturally suited to carry out predictions in message-passing neural network schemes. In bio-/chemo-/materials informatics, the challenge is to identify or design algorithms capable of obtaining molecular properties from input graphs and to follow these properties in time. We have developed a series of 2D-graphs, at various levels of granularity of re\-presentation, and associated algorithms in order to analyze physical and chemical structures and properties from atomistic molecular dynamics (MD) simulations (DFT-based MD and classical force field FF-MD). In these developments, our aim was to ensure that 2D-graphs and algorithms could be applied without any modification to ``simple'' isolated molecules, as well as to assemblies of molecules and to more complex liquid and solid states of matter, including inhomogeneous solid/liquid interfaces~\cite{bougueroua2018graph,bougueroua2023algorithmic,MolPhys2023,JCTC_ylene_2024,hashemi2022renegate,Sana-JMolStruct-2018,hashemi2023hirex}. This paper reviews some of our developments and achievements, which are also included in the GaTewAY software~\cite{bougueroua2022gateway,MolPhys2023,bougueroua2023algorithmic}. Other research groups, also experts in MD simulations, have worked on 2D-graphs. In the last decade, there have been developments of algorithmic graph theory devoted to the analysis of various types of molecular dynamics simulations, from, e.g., the conformational analysis of gas phase molecules and clusters, to their chemical reactivity, to the dynamics of H-bonds in liquids, to the dynamics of the solvation shells of ions in liquids, to the structural and dynamical analysis of complex aqueous interfaces in condensed matter~\cite{mooney_molecularnetworks_2012,ozkanlar_chemnetworks_2014,han2018graph,tenney2013analysis,pietrucci2011graph,hudelson_determining_2012,pietrucci2014fate,martinez2015automatedb,martinez2015automated,Vasudevan_JPCB_2021,bougueroua2018graph,Sana-JMolStruct-2018}. Section~\ref{sec:meth} of this paper reviews the definitions used in the atomistic 2D-MolGraphs developed by our group~\cite{bougueroua2018graph,Sana-JMolStruct-2018}, and the associated algorithms, including the key isomorphism algorithm ensuring that a whole trajectory can be analyzed in terms of the relationships between the conformations explored over time. The topological analyses hence provide a statistical view over the whole timescale of the trajectory and over the whole set of conformations explored. Statistics is crucial for the detailed knowledge of the dynamics of isolated molecules as well as for the\unskip\break dynamics of molecules in the liquid state, as applications in Section~\ref{sec:apps} will show. We also review our recent developments of 2D-graphs consisting of H-bonded cycles and the key polymorphism algorithm~\cite{JCTC_ylene_2024}, which have been built to go beyond the atomistic representation in 2D-MolGraphs and to be able to represent \mbox{(bio-)molecules} whose 3D structures and dynamics are solely based on hydrogen bonds. The algorithms we have developed enable us to track the complex conformational dynamics of flexible H-bonded molecules in real time. Section~\ref{sec:apps} will show that, while the interpretation of the complex conformational dynamics of a highly flexible hexapeptide in the gas phase would remain elusive at the atomistic level of representation (2D-MolGraphs), it is well understood by means of coarse-grained graphs of H-bonded cycles (polygraphs) and by means of polymorphic ``metastructures''. Only the coarse-grained representation in the graphs allows such comprehension. Section~\ref{sec:apps} will further show how atomistic 2D-MolGraphs can easily be included in automated high-throughput in silico reactivity workflows and how essential they are in some of the decisive steps to be taken in these workflows. We implemented the 2D-MolGraphs in the computational catalytic reaction space exploration method ReNeGate~\cite{hashemi2022renegate} and the high-throughput reactivity screening HiREX workflow~\cite{hashemi2023hirex}, specifically designed to explore realistic catalytic systems and identify thermodynamically feasible chemical transformations, corresponding to secondary catalyst deactivation and inhibition paths. Prospects and new developments in progress are discussed in Section~\ref{sec:concl}, which, with this review of methods and applications, we hope will spark further requests for additional algorithms and new applications in our physical chemistry and chemistry communities. \section{Methods}\label{sec:meth} \subsection{2D-MolGraph for modeling a conformation}\label{sec:2Dmolgraph} Our developments have been aimed at defining topological 2D molecular graphs (labeled 2D-MolGraph) and associated algorithms in order to automatically analyze MD trajectories from the knowledge of the time evolution of the conformations and hence automatically detect conformational changes through topological changes. Our 2D-MolGraphs share a similar degree of granularity in representing the topology of molecular systems to the one used by previous implementations in the literature in chemistry~\cite{ozkanlar_chemnetworks_2014,mooney_molecularnetworks_2012,pietrucci2011graph}, i.e., a vertex in the 2D-MolGraph represents an atom or a molecule and an edge between two vertices represents the interactions/bonds (covalent bond, hydrogen bond\,\ldots) between two atoms/molecules. Most of the literature on graphs in the chemical community does not consider the chemical nature of the atoms in the vertices, which is not efficient for recognizing identical structures where chemically identical atoms have been swapped. Moreover, they lack specific chemical information (e.g., covalent bonds, hydrogen bonds, exchange of atoms in homogeneous clusters, etc.) that might be relevant for a more detailed characterization of the structures. One crucial step in our method has been to define a model that represents any molecular conformation with the right level of granularity and be transferable without any modification from gas phase molecules and clusters to the condensed phase (solids, liquids, interfaces between solids and liquids). To that end, we have chosen to define any molecular conformation by a colored mixed graph $G=(V,E_{\mathrm{C}},A_{\mathrm{H}},E_{\mathrm{I}},E_{\mathrm{O}} )$, with both (directed) arcs and (undirected) edges. Such a graph is denoted 2D-MolGraph, in which each vertex represents an atom while the edges represent covalent bonds ($E_{\mathrm{C}}$), hydrogen bonds ($A_{\mathrm{H}}$), ionic (or electrostatic) interactions ($E_{\mathrm{I}}$) typically between a cation/anion atom and other atoms, organometallic interactions ($E_{\mathrm{O}}$) between metallic atoms and their surrounding. Only the hydrogen bonds are associated to directed edges (from the donor to the acceptor atom). The hydrogen atoms in a molecular system are not included in vertices of the 2D-MolGraph. Instead, their presence is solely known by directed edges. Hence, any hydrogen that is not involved in a hydrogen bond is not represented in the 2D-MolGraph. \begin{figure*} \includegraphics{fig01} \vspace*{3pt} \caption{\label{fig:grmixt}Illustration of the passage from a 3D conformation (a) to a 2D-MolGraph (b). (a)~One snapshot in the 3D representation of the gas phase Z-Ala$_6$--COOH peptide ($\mathrm{C}_{26}\mathrm{H}_{39}\mathrm{N}_{7}\mathrm{O}_{8}$) extracted from a MD trajectory. Carbon atoms are colored turquoise, nitrogen atoms blue, oxygen atoms red, and hydrogen atoms light gray. (b)~Associated topological 2D-MolGraph. Vertices are colored dark gray, blue, and red, corresponding to carbon, nitrogen and oxygen atoms, respectively. Hydrogen atoms are not included in vertices, their knowledge is included only through directed edges that represent hydrogen bonds. Edges are black lines for a covalent bond and red dashed lines (arcs) for hydrogen bonds, the latter directed from the donor to the acceptor of the H-bond.} \vspace*{5pt} \end{figure*} The definition of bonds and interactions is mainly based on Euclidian distances. The Euclidian distance between a pair of atoms $[a, b]$ with respective Cartesian coordinates $(x_a, y_a, z_a)$ and $(x_b, y_b, z_b)$ is: $\sqrt{(x_a - x_b)^{2}+(y_a - y_b)^{2}+(z_a - z_b)^{2}} $. There is a covalent bond or an interaction between two atoms if the Euclidean distance is less than a cutoff distance $D_r$. For covalent bonds, the algorithm defines the $D_r$ distance by the sum of covalent radii of atoms $a$ and $b$ with an additional 2\% margin {(that typically includes the effect on distances from vibrational motions)}. For hydrogen bonds, the algorithm sets up the default $D_r$ value between the hydrogen atom (donor) and the acceptor atom (heavy atom) to $2.3$~\AA, which can be changed by the user. For the organometallic and ionic interactions, the user is free to set case-specific $D_r$ distances. For example, the distance between manganese and oxygen atoms used in one of the applications in Section~\ref{sec:apps} was set to $2.44$~\AA. This choice was made because the developers assume that the covalent bonds are stronger than the other types of interactions between atoms. More details are found in~\cite{bougueroua2018graph,bougueroua2022gateway,hashemi2022renegate}. \begin{figure*} \vspace*{-1pt} \includegraphics{fig02} \vspace*{-1pt} \caption{\label{fig:matrix-adjacence}Illustration of an adjacency matrix for a selected portion of the Z-Ala$_6$--COOH peptide shown in Figure~\ref{fig:grmixt}a (3D-structure) and Figure~\ref{fig:grmixt}b (2D-MolGraph).} \vspace*{-1pt} \end{figure*} One can easily define and implement new relevant interactions that are needed to describe a given molecular system, and hence augment the number of definitions for the edges in the 2D-MolGraphs. In order to take into account the chemical type of the atoms in a 2D-MolGraph, we apply a special case of graph coloring, such that the vertices of a given 2D-MolGraph display the same color \textit{if and only if} the corresponding atoms have the same chemical type (see Figure~\ref{fig:grmixt}). Figure~\ref{fig:matrix-adjacence} illustrates an adjacency matrix built prior to the construction of a 2D-MolGraph. The matrix shown here is associated to a selected part of the Z-Ala$_6$--COOH peptide from Figure~\ref{fig:grmixt}a (3D structure) and Figure~\ref{fig:grmixt}b (2D-MolGraph). As the peptide contains 80~atoms, only a selected part of the peptide has been extracted here for this illustration. The figure shows that covalent bonds and hydrogen bonds of interest in the conformation of this peptide are encoded in the adjacency matrix (with crosses and circles, respectively) and that the matrix has the colored information of the actual chemical nature of the atoms. The graph on the right side of the figure is the subgraph of the 2D-MolGraph associated to the selected part of the molecule encoded in the matrix. With 2D-MolGraphs in our hands, the exploration with time of molecular conformations along MD simulations can easily be seen as the exploration of graph topologies, that can be tracked using graph theory-based methods, such as isomorphism. Graph isomorphism as defined in~\cite{mckay1981practical} allows representation of each conformer with a fingerprint graph and comparisons between graphs. In our applications, isomorphism consists in comparing the distribution of edges between two 2D-MolGraphs: if the graphs compared have the same set of bonds/interactions connected to the same set of atoms (in terms of chemical types/colors for the graphs), these graphs are then isomorphic, i.e., the two graphs are identical. More formally, an isomorphism is a bijection between the vertex sets of the two graphs if and only if it induces a corresponding bijection between their edge sets (if such an isomorphism exists, the two graphs are said to be isomorphic). A (non-polynomial) algorithm to check if such an isomorphism exists is proposed in~\cite{luks1980isomorphism,mckay1981practical,McKay201494,hartke2009mckay}. Isomorphism checks, together with keeping the chemical nature of the atoms, are the key components of the conformational search over MD trajectories. The changes in conformations are hence followed over time by scanning trajectories for changes in bonding patterns of choice (among hydrogen bonds, proton transfers, coordination numbers, covalent bonds, and organometallic interactions). Once the different conformations of the molecular system have been found by graph analysis, a graph of transitions can be generated, similar to the ones considered in different analysis or generation of temporal graph sequences~\cite{Casta-Dag-21,Liu-KDD-23,Plmaper-EST-23}. This graph has its vertices composed of the conformations that have been identified and its edges are composed by the transitions found between conformations. Both vertices and edges in the graph of transitions contain information on the percentage of existence of a conformation over the duration of the trajectory (for the vertices) and the percentage of times a transition between two vertices/conformations has been seen (for the edges). In one glance, one can hence see the relationships between the conformations and associated statistics. \begin{figure*} \includegraphics{fig03} \caption{\label{fig:grphtrans}Graph of transitions. See text for nomenclature and colors.} \end{figure*} Figure~\ref{fig:grphtrans} shows an example of a graph of transitions. This graph is composed of four vertices, each vertex represents one molecular conformation that has been identified by isomorphism along the trajectory. Each vertex has a number and a label that represents the conformation. For instance, one vertex in this picture is labeled with ``1'' and ``N1-O1'', meaning this is the first conformation identified over the trajectory, and the associated conformation has one hydrogen bond between atoms ``N1'' and ``O1''. The ``68.97\%'' number in the vertex is the percentage of time this conformation has been seen over the whole trajectory. Hence, in this graph of transition, conformation 1 is the most frequent conformation that has been observed over the trajectory, with a total percentage of appearance of ${\sim}$69\% over the trajectory. The edges between the vertices are labeled with two kinds of information: (1)~the total frequency rate for going from one conformation to the other one, (2) the bond(s)/interactions(s) that have changed when going from one conformation to the other. For instance, one can observe a large conformational dynamics between conformations $1$ and $2$ in the graph of transitions in Figure~\ref{fig:grphtrans}; such an event occurs around $130$ times, back and forth. Also of note, the hydrogen bond N1$\cdots$O1 disappears when going from conformation~1 to conformation~2 (and appears on the reverse way). Conformation~2 in the graph of transitions is labeled with ``2'' and ``.'', the latter meaning that this conformation contains no hydrogen bond. Conformation~2 has been seen ${\sim}$30\% over the whole trajectory. The colors of the vertices in the graph of transitions directly give the most relevant conformations in terms of appearance periods. We hence colored red the conformations that appear at least $P_{\mathrm{min}}$\% (an input parameter that the user can change the default value of 4\% has been used here) and the ones in green occur below this threshold. All the conformations explored along the MD simulations can be kept in the graph of transitions. Such information might indeed be useful for some analyses, typically when rare events (rare conformations) are investigated. The user can modify this at will. \subsection{From a topological 2D-MolGraph to a coarse-grained graph of H-bonded cycles}\label{polygraphs} In~\cite{MolPhys2023,bougueroua2023algorithmic} we have shown that the rationalization of the conformational dynamics becomes complex and almost impossible to achieve for flexible H-bonded molecules that isomerize frequently over time through the dynamics of breaking and forming of their H-bonds. As shown in these references, a simple short peptide such as \zala\ (illustrated by one \molgraph\ in Figure~\ref{fig:grmixt}) already shows a high flexibility of its network of N--H$\cdots$O H-bonds at relatively low temperatures (gas phase MD trajectories). Numerous breaking/forming of H-bonds were observed over the trajectories, which signaled the appearance and disappearance of several conformers of the gas phase peptide numerous times along the trajectory. This high dynamical flexibility of the H-bond network, however well represented by the 2D-MolGraphs and their graph of transitions~\cite{MolPhys2023,bougueroua2023algorithmic}, prevents a clear rationalization of the actual conformational dynamics of the peptide. In particular, some of the H-bonds in the network, though formed between different atoms, seem to play a similar role into the final 3D structure of the peptide molecule, i.e., into the final folded/semi-folded/unfolded skeleton. This is not captured by the 2D-MolGraphs at the atomic level of representation, but this similarity of H-bonded cycles can be captured by graphs defined at a higher/coarser granularity of representation. We therefore defined graphs in which the vertices are directly associated to the H-bonded cycles formed, whose polymorphism is furthermore taken into account. These concepts are now explained. Cycles in molecular graphs have been shown to be good representations of the structure of molecular systems~\cite{Berger2017,GLEISS2001,Ilemo2019}. We hence proposed a representation of each conformation identified by a 2D-MolGraph over a trajectory based on a well chosen set of cycles. The cycles of interest to us are the ones formed by at least one H-bond. With these, we will be able to quantify and follow in time the changes in the H-bonded network of molecular structures. Our current developments have been done for gas phase molecules only~\cite{JCTC_ylene_2024}. Note that the number of cycles in a graph can be exponential with respect to the number of vertices. Therefore in our approach, we only consider a subset of cycles in the \molgraph, called ``minimum cycle basis'', restricted to the H-bonded cycles for each \molgraph\ of the trajectory. Given a graph, finding a minimum cycle basis, which is not necessarily unique, can be done in polynomial time with an evolution of the algorithm of Horton~\cite{Horton1987}. More details can be found in~\cite{JCTC_ylene_2024}. Given a \molgraph, we compute its minimum basis set of H-bonded cycles. The associated graph of cycles is defined as a graph in which the set of vertices is the cycle basis, and there is an edge between two vertices (cycles) if and only if these two cycles interact, i.e., they share \textit{at least} one covalent bond or one hydrogen bond in the \molgraph. Hence each 3D-conformer in the trajectory of a molecule is\unskip\break represented by a graph of cycles. The evolution in time of the conformations in a MD trajectory can be represented by the sequence of their graphs of cycles. Figure~\ref{fig:step1_conf} illustrates the transformation from a topological graph (Figure~\ref{fig:step1_conf}a) to a graph of cycles (Figure~\ref{fig:step1_conf}c). Considering the \molgraph\ in Figure~\ref{fig:step1_conf}a and the minimum cycle basis in Figure~\ref{fig:step1_conf}b, there are four H-bonded cycles, each of them of various size, composed of one or several hydrogen bonds. Hence, the pink vertex in Figure~\ref{fig:step1_conf}c is built on one H-bonded cycle composed of six vertices/atoms in the \molgraph\ (this is a six-membered H-bonded cycle), while the larger orange vertex/cycle is built upon two hydrogen bonds (see the two directed dashed red edges in Figure~\ref{fig:step1_conf}a). Figure~\ref{fig:step1_conf}c shows the graph of cycles obtained from the minimum cycle basis shown in Figure~\ref{fig:step1_conf}b. In this graph, the vertices are the H-bonded cycles, labeled by the heavy atoms involved in the hydrogen bond(s) producing them. Taking once again the examples of the pink and orange cycles/vertices depicted in Figures~\ref{fig:step1_conf}b--\ref{fig:step1_conf}c, the pink vertex is labeled N3O3 as it is built on the N3--H$\cdots$O3 hydrogen bond, while the orange vertex is labeled with its two constitutive H-bonds N1--H$\cdots$O8 and N5--H$\cdots$O4. The orange, pink and blue vertices/cycles interact with each other as seen through the edges connecting these three vertices (i.e., sharing at least one covalent bond or one hydrogen bond). On the other hand, the green vertex (related to the H-bonded cycle N7O7) interacts only with the orange one. \begin{figure*} \includegraphics{fig04} \caption{\label{fig:step1_conf}Example for the coarse-grained representation of a \molgraph\ for one conformation of the gas phase peptide \zala. See text for details.} \end{figure*} Given a MD trajectory, each \molgraph\ is now associated to a minimum cycle basis and to the corresponding graph of cycles. In other words, there is one graph of cycles per 2D-MolGraph (i.e., per identified molecular conformation). The set of conformational isomers explored over the MD trajectory can furthermore be summarized by one single graph of cycles uniting those of all the identified conformers. This union takes into account all the possible H-bonded cycles as well as all the possible interactions between these cycles which were observed in all the identified conformers. Some of these H-bonded cycles can be identical in different conformers. Some of these H-bonded cycles can be similar to each others in the sense that they are built upon different donor/acceptor atoms but they are playing the same role in the final structure of the molecule. One therefore has to recognize the similarity between\unskip\break the H-bonded cycles of the different miminum cycle bases (i.e., similarity between the H-bonded cycles in the different identified conformations). To that end, the next stage was to develop an algorithm that can group together the cycles of different minimum cycle bases which are similar, i.e., the ones playing the same role in the structure of different conformers. The similarity constraint that we introduced is the basis to cluster the union of cycles in all the minimum cycle bases of all the conformers of the trajectory that we finally want to obtain. Each part in such a clustering is called a \textit{polymorphic cycle}, in which all cycles are considered as different forms of a same cycle in the structures of the conformers it appears in. The hypothesis that we made is that a cycle can evolve over time, that is to say that its set of links can evolve while retaining its same role in the molecular structure. The trajectory can therefore be seen as the interaction of the H-bonded cycles evolving over time in their atomic structure (polymorphism) and still interacting in the same way, but with some of these cycles appearing or disappearing over time. The final graph thus obtained is called a polygraph, a contraction of ``polymorphic cycle graph''; note that this definition is different from the one of a polygraph as a generalization of the directed graph~\cite{Babic-DM-86}, or that used in chemoinformatics related to polymers~\cite{Juvana-DM-97}. The whole algorithms have been detailed in~\cite{JCTC_ylene_2024}. We now show an application of this methodology for the gas phase \xchondro\ molecule, for which a 3D snapshot is reported in Figure~\ref{fig:CS2S4S_3D}. Figure~\ref{fig:cs2s4sPM6_polygraphFull} reports the polygraph obtained for a 600~K trajectory of this gas phase molecule. Each vertex of the polygraph is labeled by ``$\mathbf{P}$-XX'' where $\mathbf{P}$ stands for polycycle and XX is the set of oxygen and/or nitrogen atoms (with their associated label number in the list of the molecule's atoms) involved in the H-bond of the polycycle (i.e., one of the two partners which close the H-bonded cycle). Each polycycle has been assigned a given color. Note that a polycycle can be built over one or over several H-bonds (two to three for this molecule). In that case, the atoms involved in the series of H-bonds are written on a line with a ``,'' between them. See some examples in the vertex colored in violet for several instances of multiple H-bonds forming a polycycle. \begin{figure} \includegraphics{fig05} \caption{\label{fig:CS2S4S_3D}A 3D representation of \xchondro\ (C$_{14}$H$_{20}$NO$_{18}$S$_{2}$). Colors for the atoms: carbon in turquoise, nitrogen in blue, oxygen in red, hydrogen in light gray, and sulfur in yellow.} \end{figure} \begin{figure} \includegraphics{fig06} \caption{\label{fig:cs2s4sPM6_polygraphFull}The polygraph obtained for a MD trajectory of the gas phase \xchondro\ molecule.} \end{figure} Figure~\ref{fig:cs2s4sPM6_polygraphFull} shows that this polygraph is a complete graph made of five vertices. It is complete because there is an edge between each pair of vertices. As a reminder, there is an edge between two polycycles whenever these two cycles have atoms that share at least one covalent bond or one H-bond in the initial 2D-MolGraph (i.e., these polycycles are interacting with each other within the molecular structure). Because the polygraph is complete, the number of polycycles cannot be reduced. All vertices involve different polymorphic identities. Four of the five vertices have six or seven different actual identities for the H-bonded cycle that forms the polycycle (dark green, pink, green, and orange), while there is an even larger diversity of 14 different identities for polycycle $\mathbf{P}$-O17 (purple). For instance, vertex $\mathbf{P}$-O7 is associated to a polycycle that can have six different identities of the H-bond that closes the cycle/polygon with the O7 atom. As can be seen, this O7 can either H-bond to the H atom carried by the O13, or to O9, O16, O7, O15 and O14. In the vertex in the violet color, some of the lines report other identities than the $\mathbf{P}$-O17 identity. For instance, one line reports O1-O2 and O2-O17, which means that the $\mathbf{P}$-O17 polycycle is built over two simultaneous hydrogen bonds, i.e. the O2-O17 H-bond but also the O1-O2. \section{Review of selected applications}\label{sec:apps} \subsection{H-bond dynamics of a flexible gas phase peptide and the relevance of the coarse-grained graphs of H-bonded cycles} The conformational dynamics of the gas phase hexapeptide \zala\ (C$_{26}$H$_{39}$N$_7$O$_8$, $80$ atoms, 3D illustration in Figure~\ref{fig:Z-Ala6_3D}) is analyzed hereby both in terms of the 2D-MolGraphs (atomistic level of topology representation) and in terms of the coarse-grained polymorphic cycles (coarse-grained level of topology representation) in order to present the strengths and limitations of each level. An ab~initio MD (AIMD) trajectory of ${\sim}$6.0~ps (12,294 snapshots, $\updelta t$ of 0.5~fs) at ${\sim}$450~K serves as the basis for the conformational dynamics to analyze. This temperature has been chosen as a good example for the analysis of trajectories where several medium and large conformational changes are expected to occur, which are always rather hard to characterize without the help of topological graphs. \begin{figure} \includegraphics{fig07} \caption{\label{fig:Z-Ala6_3D}A 3D representation of the hexapeptide \zala. Colors for the atoms: dark gray for carbon, dark blue for nitrogen, red for oxygen, white for hydrogen.} \end{figure} \subsubsection{A high conformational dynamics of \zala\ provided by the 2D-MolGraphs} The analysis of the trajectory in terms of topological molecular 2D-MolGraphs~\cite{bougueroua2018graph,MolPhys2023} indicates that $93$ different conformations are sampled over ${\sim}$6~ps (which hence shows a rather highly dynamical peptide), built over nine different H-bonds. The structure of each identified conformer is composed of one to six H-bonds that are formed simultaneously. The peptide is found in either opened structures where a low number of simultaneous H-bonds are present or in H-bonded folded structures. As illustrations, Figure~\ref{fig:confabzala} presents two 2D-MolGraphs corresponding to two identified conformers of \zala, respectively produced by five (Figure~\ref{fig:confabzala}c) and three H-bonds (Figure~\ref{fig:confabzala}d). Two of these H-bonds (N7$\cdots$O7 and N6$\cdots$O5) are present in both conformations. \begin{figure*} \includegraphics{fig08} \vspace*{-1pt} \caption{\label{fig:confabzala}Two 2D-MolGraphs (c and d) of the gas phase peptide \zala\ and their associated 3D structures (a and b) extracted from the MD trajectory analysis. Colors of the vertices: dark gray for C atoms, dark blue for N, red for O. No vertices for hydrogens, their knowledge is in the directed edges (arcs). Solid black edges in the graphs are for covalent bonds, the red arcs are the hydrogen bonds directed from the donor to the acceptor (heavy) atoms.} \vspace*{-3pt} \end{figure*} Figure~\ref{fig:graphtranszala} is the graph of transitions that summarizes the whole 6~ps of conformational dynamics of \zala\ at 450~K. As can be immediately seen, this graph of transitions is composed of an extremely high number of vertices and edges connecting these vertices, which is the signature of the high dynamicity and high flexibility of the peptide. There are $93$ vertices in the graph for the $93$ different conformers of \zala\ found. In the graph of transitions reported in Figure~\ref{fig:graphtranszala}, most of the identified conformers (in the green vertices) appear over very short periods of time (less than 0.24~ps) while two conformers (red vertices) appear over larger durations\unskip\break of time. \begin{figure*} \includegraphics{fig09} \caption{\label{fig:graphtranszala}Graph of transitions for the 450~K dynamics of the gas phase \zala\ peptide. The vertices indicate the explored conformers. Vertices in red are for conformers with a total percentage of appearance $P_{\mathrm{min}}$\% greater than 4\% of the dynamics time (parameter $P_{\mathrm{min}}$\% can be modulated at will), vertices in green are for conformers with $P_{\mathrm{min}} \% < 4\%$. Directed edges between vertices indicate transitions between two conformers as observed over time. The labels on each edge provide the total percentage of occurrence of the transition and the associated chemical change(s) that occur. A zoom over a small portion of the graph of transitions is provided on the left-hand side.} \end{figure*} While the high number of vertices and edges in the graph of transitions in Figure~\ref{fig:graphtranszala} illustrates the high flexibility of the peptide at 450~K that is nicely\unskip\break captured by our 2D-MolGraph topological graphs, it also illustrates the limit of the atomic level of granularity used in the 2D-MolGraphs for highly flexible molecular systems where presumably too much information is provided in the graph of transitions and is hard to process. Furthermore, some of the periods of times that the 2D-MolGraphs analyze as breaking/forming covalent bonds and/or H-bonds correspond in practice to the dynamics of these bonds around the threshold values employed in the method for conformational recognition. Some of the transitions observed between conformers are thus in practice the actual signature of the existence of one single ``meta-conformation'' around which dynamicity and flexibility occur. This information can however not be extracted from the atomistic topological 2D-MolGraphs. To get that information, one has to analyze the trajectory with a coarse-grained topology representation, whose results are described in the following section. \subsubsection{Polygraphs are the good coarse-grained representation to analyze the conformational dynamics of \zala} To go beyond the limitations of the atomistic representation in the topology highlighted above in the case of the highly flexible \zala, we now apply the coarse-grained topology representation described in Section~\ref{polygraphs} to analyze the conformational dynamics of this peptide: it now consists in representing a molecule through the ensemble of its H-bonded cycles and to apply an algorithm of polymorphism in order to recognize the cycles that are isomorphic to each others. Figures~\ref{fig:zala6450DFT_polygraph} and~\ref{fig:zala6450DFT_chrono} illustrate the two immediate outputs of this analysis, respectively showing the global polygraph generated over the whole 450~K trajectory in Figure~\ref{fig:zala6450DFT_polygraph} and the chronogram (time evolution) of the nine identified polymorphic cycles in Figure~\ref{fig:zala6450DFT_chrono}. \begin{figure*} \includegraphics{fig10} \vspace*{-4pt} \caption{\label{fig:zala6450DFT_polygraph}The global polygraph obtained over the 6~ps AIMD trajectory of the \zala\ peptide at 450~K. The color coding of the vertices is identical to the color coding of the lines in the chronogram in Figure~\ref{fig:zala6450DFT_chrono}. Each vertex of the polygraph is labeled by ``{P}-XX'' where {P} stands for polycycle and XX is the set of the oxygen atom and/or the nitrogen atom (with their associated label number in the list of the atoms of the molecule) involved in the H-bond of the polycycle. Each vertex contains the list of the polymorphic cycles, each one being labeled by the atoms' names and labels in the H-bond(s) building the H-bonded cycle.} \vspace*{-4pt} \end{figure*} There is a total of $23$ H-bonded cycles over the $93$ 2D-MolGraphs extracted from the ${\sim}$6~ps trajectory of \zala\ at 450~K. Once polymorphism has been applied over these $23$ H-bonded cycles, some of the cycles could be merged, thus resulting into nine polymorphic cycles. Each polycycle contains between one and nine possible identities. The obtained global polygraph that represents the whole 6~ps trajectory is shown in Figure~\ref{fig:zala6450DFT_polygraph}, built over the nine vertices that represent the nine identified\unskip\break polymorphic H-bonded cycles. One can therefore see that five of the vertices/polycyles have only one possible identity in terms of the H-bonded cycle (note that the H-bonded cycle for the polycyle P-N1,O8,N2 is built upon two hydrogen bonds, i.e., N1$\cdots$O8 and N2$\cdots$O8), two other vertices have a limited number of isomorphic identities (two identities for the dark blue vertex P-O4,N5, three identities for the magenta vertex P-N3). Two vertices are much more polymorphic as they adopt four (cyan vertex P-O4,N5) and up to nine (orange vertex P-O8) different polymorphic identities. Furthermore, one can remark that several of these polycycles are built upon two or three H-bonds. For instance, in the orange P-O8 polycycle, the first identity ``N1-O8,O4-N4'' is built upon two H-bonds (N1$\cdots$O8 and O4$\cdots$N4), the first identity of the cyan polycycle P-O4,N5 is built upon three H-bonds (O2$\cdots$N3, N1$\cdots$O8, and O4$\cdots$N5). \begin{figure*} \includegraphics{fig11} %\vspace*{-4pt} \caption{\label{fig:zala6450DFT_chrono}Chronogram obtained over the 450~K AIMD trajectory of the \zala\ peptide. Five time periods have been identified (see text for details). The color coding of the lines is identical to the color coding of the vertices in the polygraph in Figure~\ref{fig:zala6450DFT_polygraph}.} \vspace*{2pt} \end{figure*} The low/high number of identities within each polymorphic H-bonded cycle informs on the sections of the peptide with a low/high structure\unskip\break flexibility. The higher polymorphic nature of two of the polycycles in \zala\ corresponds to high structure flexibility in these two zones, while maintaining the same general role of the H-bonded cycle into the final 3D structure. This high flexibility was already observed in the previous section with the atomistic topological analysis of the 2D-MolGraphs. The analysis with the polymorphic H-bonded cycles immediately identifies the underlying ``meta-structures'' that are of interest for the comprehension of the conformational dynamics of the peptide. \looseness=-1 In the polygraph, there is an edge between two vertices whenever two polycycles interact together, i.e., the conformations belonging to each vertex share at least one covalent bond or one H-bond. The polygraph in Figure~\ref{fig:zala6450DFT_polygraph} is not complete, i.e., not all the polycycles are directly connected in pairs by an edge. This incompleteness can be explained by the following two reasons that are related to the rules applied for building a polygraph: (i)~either there is a\unskip\break conformation for which the identities of the two polycycles appear simultaneously, without interaction (the criterion for interaction/edge is the sharing of at least one covalent bond or one H-bond); or (ii)~there is no conformation for which the two identities appear simultaneously, however their merging does not respect the polycycle rules. For example, we found that the blue P-O4,N5 polycycle/vertex and the orange P-O8 polycycle/vertex do not interact because there is at least one identity from P-O4,N5 and one identity from P-O8 that appear simultaneously in one conformation of \zala. As these two polycycles do not share one covalent bond/H-bond, there is thus no interaction/edge in the polygraph. The chronogram of the \zala\ polycycles, presented in Figure~\ref{fig:zala6450DFT_chrono}, shows the evolution in time of the nine polycycles at the temperature of 450~K. The reader has to keep in mind that each polycycle can have several identities in terms of the atoms that form the H-bonded cycle (see discussion above). Five distinct periods of time can be extracted, labeled as period~1 to period~5 in the figure. Each time period is associated to a different structural organisation of the nine polycycles in \zala. During period 1, six polycyles coexist. Period~2 starts as the polycycle P-N3 disappears and polycycle P-O2,N3 appears. There is then the simultaneous disappearance of this polycycle at time ${\sim}$2.2~ps and appearance of polycycle P-N1,O8,N2 at the slightly earlier time ${\sim}$1.9~ps, which marks a clear change in the 3D structure of the peptide, and thus the end of period~2 and start of period~3. Interestingly, polycycles P-O7,N7 (green line), P-O4,N5 (dark blue) and P-O4,N5 (cyan) are formed over the whole 6~ps trajectory (or almost always, with smallish periods of time of interruption). These H-bonded polycycles thus form strong pillars for the global 3D structure. See more details in~\cite{JCTC_ylene_2024}. One has to remark that one of these three polycycles (P-O4,N5, darker blue) can adopt up to four different identities. This is an important pillar of the 3D structure with high flexibility. Beyond the time evolution, one further crucial information extracted from the chronogram is which polycycles can be formed simultaneously and which ones cannot be present simultaneously in the 3D conformation of the peptide. For example, it is possible to form the polycycle P-O4,N4 (light blue) at the start of period~5 only if none of the P-O4,N5 (cyan) and P-O5,N6 (purple) polycycles are present. P-O4,N4 seems to coexist only with P-O2,N3 (light pink) and P-N1,O8,N2 (dark orange). Polycycles P-O4,N4, P-O2,N3, P-N1,O8,N2, P-O5,N6, and P-N3 need hence certain structural conditions for them to be exist. This can be explained by the location of H-bonded polycyle P-O4,N4 in regard to polycycles P-O4,N5 (both of them) and P-O5,N6 within the 3D structure of \zala. The polygraph made of nine vertices and their pairwise connected edges in Figure~\ref{fig:zala6450DFT_polygraph} gives the global/statistical view of the whole 6~ps trajectory in terms of a ``global metastructure'' of the peptide. However, the actual details of the dynamics are lacking in this global polygraph. This could be inferred from the chronogram in Figure~\ref{fig:zala6450DFT_chrono} where we already saw that not all the polycycles/vertices of the global polygraph could coexist simultaneously over time. \begin{figure*} \includegraphics{fig12} %\vspace*{-4pt} \caption{\label{fig:z6dft1}Sequence in time of the sub-polygraphs of the H-bonded polycycles of \zala. See text for all details and comments of this chart.} \vspace*{-5pt} \end{figure*} Figure~\ref{fig:z6dft1} now presents one illustration of the evolution with time of the sub-polygraphs from the global polygraph over a short time period, i.e., between times 2.19~ps and 2.51~ps (period~2 and period~3 in Figure~\ref{fig:zala6450DFT_chrono}). The sequence of sub-polygraphs is reported with the following conventions of colors: the polycycles that are present among the nine\unskip\break possible polycycles are colored in green while the ones that are absent are colored in dark gray. The arrows between the sub-polygraphs can be colored in black, in blue, or in red. The two arrows in red mark the sequence of starting time of the period under scrutiny (2.19~ps) and ending time (2.51~ps), respectively. The arrow in black marks ``simple transitions'' between sub-polygraphs/3D-structures/conformations of \zala. These arrows can be found forward (${\rightarrow}$) and backward (${\leftarrow}$), as there are multiple forward/backward isomerisations over time. Marking each one of these transitions over the ``real time'' would be too messy and would make the plot too cumbersome. The blue arrow marks a distinct transition in the sequence of sub-polygraphs, i.e., in the time sequence of the 3D conformations adopted by the peptide over this period of time. There is indeed no transition observed back to structures adopted previously (in time). Therefore, we marked the whole chart by two distinct rectangles, one outlined with red dots (top of the figure) and one with gray dots (bottom). In the gray rectangle, one can now see the sequence between the sub-polygraphs by following the dark arrows, some of them going in one direction only, others with back/forth directions. Some of the key moments in time are indicated in blue over certain dark arrows. \begin{figure*} \includegraphics{fig13} \vspace*{4pt} \caption{\label{fig:aw_gr_br_conf}Illustration of one 2D-MolGraph per hydrophobic aqueous interface: air/liquid water ({left}), graphene/liquid water ({middle}), BN/liquid water ({right}). Only the water molecules located in the BIL~\cite{Pezzotti_PCCP_2018,pezzotti2021molecular} are taken into account for the 2D-MolGraph analysis. Vertices of the graph represent the oxygen atoms of the water molecules (red); the dashed red arcs represent the H-bonds between two water molecules oriented from donor to acceptor.} \vspace*{6pt} \end{figure*} The conformation dynamics seen in this figure over a very short period of time nicely illustrates the complex dynamics in the H-bonds of \zala\ at 450~K, and the high flexibility of the peptide that we have already discussed several times above. Figure~\ref{fig:z6dft1} also nicely shows the fast exchange between polycycles. \subsection{Topological 2D-MolGraphs easily capture complex condensed phase H-bonded networks} Condensed phase, where liquid water is present, is another area where the 2D-topological graphs can be of great help in understanding H-bond networks. Here we present the use of the atomistic 2D-MolGraphs in order to unravel the structure of water in inhomogeneous molecular systems made of an interface between liquid water and another medium. We specifically focus on three hydrophobic aqueous interfaces, i.e., air/liquid water, graphene/liquid water, and boron nitride BN/liquid water, for which we want to characterize the organization of liquid water at the interface with air or the solid. To that end, our methodology of atomistic topological 2D-MolGraphs is applied on AIMD trajectories of these three aqueous interfaces (${\sim}$50~ps time-length MD). These three interfaces have been shown to be hydrophobic by independent molecular analyses in~\cite{pezzotti2021molecular} based on a molecular descriptor of hydrophobicity developed in this latter paper and in the follow-up paper~\cite{wanlin_2023JPCL_HVdescriptor}. We have previously shown in~\cite{pezzotti20172d,Sana-JMolStruct-2018,pezzotti2021molecular} that liquid water in contact with hydrophobic surfaces forms a two-dimensional (2D) highly collective H-bonded network made by the water molecules in the layer in direct contact with the hydrophobic surface (i.e., water located in the BIL--Binding Interfacial Layer-- as defined in~\cite{Pezzotti_PCCP_2018,pezzotti2021molecular}), in which the water--water H-bonds are formed parallel to the surface. This water-collective 2D-Hbonded-Network is the molecular signature of surface hydrophobicity~\cite{pezzotti2021molecular,wanlin_2023JPCL_HVdescriptor}. Here, we illustrate the recognition of this 2D-HBonded-Network using topological 2D-MolGraphs, and how these graphs can provide details on the organization of water molecules in this collective H-bond network. Three DFT-MD trajectories have been analyzed using our graph theory algorithm: air/liquid water as the prototype of hydrophobic surfaces, graphene/liquid water, and BN/liquid water. For each trajectory, 400~snapshots were extracted and analyzed (from a total of 50~ps trajectory per system). This corresponds to roughly one snapshot every 0.1~ps of dynamics, which represents a good statistical sampling regarding the dynamics of H-bonds. The graph analyses are carried out on the BIL-interfacial region only, in which there is roughly an average of 48 water molecules (all simulation boxes are roughly equivalent in sizes). \begin{figure*} \includegraphics{fig14} \caption{\label{fig:aw_gr_br_cc}Distribution of the connected components of the 2D-MolGraphs (see text for details) for the air/liquid water interface ({left}), the graphene/liquid water ({middle}), and the BN/liquid water ({right}), obtained over ${\sim}$50~ps AIMD trajectories.} \vspace*{6pt} \end{figure*} \begin{figure*} \includegraphics{fig15} \caption{\label{fig:aw_gr_br_ring}Distribution of the size of H-bonded rings/cycles formed by the water molecules in the 2D-MolGraphs. ({left}): air/water interface; ({middle}): graphene/water interface; ({right}): BN/water interface.} \end{figure*} The 2D-MolGraphs in Figure~\ref{fig:aw_gr_br_conf} show the very specific structural property of the water molecules in the BIL of the hydrophobic interfaces: a collective arrangement of the water molecules in terms of H-bonded polygons (or rings) adjacent to each others. This collective arrangement is called the 2D-HBonded-Network~\cite{pezzotti20172d,Sana-JMolStruct-2018,pezzotti2021molecular}. To obtain the statistical view on the number of water molecules that are interconnected within the 2D-HBonded-Network, the (identified non-isomorphic) 2D-MolGraphs can be analyzed in terms of the size of the connected components, i.e., the set of subgraphs in which all vertices are connected to each others without interruption. Figure~\ref{fig:aw_gr_br_cc} illustrates the\unskip\break distribution of the connected components for the 400 2D-MolGraphs extracted for each trajectory of the three interfaces. One can immediately see that the 2D-HBonded-Network is extended over ${\sim}$90--95\% of the water molecules located in the BIL for the three interfaces, hence with a high degree of interconnectivity (very collective network between the water molecules). One can thus conclude that the water molecules in the BIL are statistically organized with the same collective HB-Network in all these hydrophobic\unskip\break interfaces. The 2D-MolGraphs provide further details on the organization of water molecules in this collective H-bonded network. One can see from the 2D-MolGraphs shown in Figure~\ref{fig:aw_gr_br_conf} that water molecules are organized in polygons/rings formed by H-bonds. Using the Horton algorithm~\cite{horton1987polynomial}, we analyzed the 2D-MolGraphs in terms of the size of H-bonded\unskip\break polygons/rings formed by the water molecules in the BIL for the three hydrophobic interfaces. The results are presented in Figure~\ref{fig:aw_gr_br_ring}. Very interestingly, though the water molecules in the three investigated BILs are assembled with the same collective 2D-Hbonded-Network, the distribution of sizes of the H-bonded polygons that build these networks are non-identical between the three hydrophobic interfaces. On the one hand, the sizes of the H-bonded polygons are centered on four to six for the air/water and BN/water interfaces. There is a clear dominant component related to H-bonded pentagons made by the water molecules at the interface with the BN surface while the formation of H-bonded tetragons and pentagons is found equivalent at the interface with air. On the other hand, the 2D-Hbonded-Network made by the water molecules at the surface of graphene is more homogeneous in terms of sizes of the polygons, where tetragons, hexagons, and heptagons have roughly the same probability of appearance, and H-bonded pentagons dominate slightly more. We hence see that water molecules predominantly form five-membered H-bonded rings/polygons at the interface with the BN surface, which can be associated to the hexagonal-templated structure of BN. The length of the C--C covalent bonds in BN is shorter than the O$\cdots$H hydrogen bonds: the best arrangement for water molecules is thus into five-membered H-bonded rings rather than six-membered rings. Moreover, we find that the percentage of water molecules giving rise to the polygons within the 2D-HBonded-Network is around 30--40\% for the three interfaces, with the following interesting ranking: one finds a larger percentage of water in the 2D polygons for the air/water interface (${\sim}$43\%) than for the graphene/water interface (${\sim}$36\%) interface, and the BN/water interface (${\sim}$34\%). Such percentages might explain the strength of the 2D-HBonded-Network found at each interface. The work in~\cite{pezzotti2021molecular} indeed showed that the strength of the 2D-HBonded-Network can be ranked as air ${>}$ graphene ${>}$ BN. In other words, the more water molecules forming rings within the 2D-Hbonded-Network (i.e., the more rings being formed), the stronger the 2D-HBonded-Network, and the more hydrophobic the interface. \subsection{Integration of 2D-MolGraphs in workflows in high-throughput in silico chemical reactivity} We further explored the utility of our 2D-MolGraph approach as the core for automated reaction network analysis workflows suitable for implementation in high-throughput in silico reactivity screening of complex multicomponent homogeneous catalytic systems~\cite{pidko2024}. We implemented the 2D-MolGraph approach in the computational catalytic reaction space exploration method ReNeGate~\cite{hashemi2022renegate} (see the workflow in Figure~\ref{fig:renegate}) and in the high-throughput reactivity screening HiREX workflow~\cite{hashemi2023hirex} specifically designed to explore realistic catalytic systems and identify thermodynamically feasible chemical transformations, corresponding to secondary catalyst deactivation and inhibition paths. \begin{figure*} \includegraphics{fig16} \caption{\label{fig:renegate}Schematic representation of the ReNeGate workflow from~\cite{hashemi2022renegate} involving the sequential reactive space exploration, the structure analysis using the 2D-MolGraphs in the fingerprinting and network construction parts, with the reaction network generation and refinement steps (trimming).} \vspace*{-2pt} \end{figure*} The method was validated by case studies on representative multicomponent (de)hydrogenation catalytic systems based on Mn(I) coordination\unskip\break complexes with a special focus on probing unconventional and less expected reaction channels, which could be responsible for loss of catalytically potent species during the initial pre-catalyst activation. An illustrative example is our computational analysis of the activation of manganese pentacarbonyl bromide (Mn(CO)$_5$Br) with inorganic alkoxide base KOiPr that is a common protocol for the experimental screening and in situ generation of Mn-based homogeneous catalysts~\cite{yang2023performance}. The reactivity exploration was carried out with parallel metadynamics simulations on a minimal model constituted by the two main reagents only (Figure~\ref{fig:MnCOBr}a) using the CREST functionality (metaD/CREST) at the GFN2-xTB level~\cite{grimme2019exploration}. The reactive trajectories were populated using the root-mean-square-deviation (RMSD) in Cartesian space as a metric for the collective variables~\cite{grimme2019exploration}, while the pushing and pulling strengths ($k$ and $\alpha$) were systematically varied over the parallel simulations. \begin{figure*} \includegraphics{fig17} {\vspace*{-2pt}} \caption{\label{fig:MnCOBr}(a)~Metadynamics simulations of the reaction between Mn(CO)5Br and KOiPr yield reactive trajectories, from which (b)~a reaction network has been constructed with the respective three representative 2D-MolGraphs presented in (c). Colors for the 2D-MolGraph vertices: dark gray for carbons, red for oxygens, dark blue for postassium, pink for the manganese, and green for the bromide. Solid black edges in the graphs represent covalent bonds, the red dashed edges represent the organometallic interactions.} {\vspace*{-2pt}} \end{figure*} The analysis of the reactive trajectories with the 2D-MolGraphs yielded a reaction network containing 12 conformers, which after trimming the edges exceeding an arbitrary threshold of 25~kcal/mol, produced the reaction network shown in Figure~\ref{fig:MnCOBr}b. State (1) corresponds to the starting configuration with unreacted Mn(CO)$_5$Br and KOiPr, which can transform to one of the new species identified by the 2D-MolGraphs from the reactive trajectories (Figure~\ref{fig:MnCOBr}c). Our automated procedure revealed that all reaction paths involve the reaction of the alkoxide nucleophile with the Mn(I)-bound carbonyl ligand to form a Mn--acyl complex (2). Due to its approximate nature, the GFN2-xTB method incorrectly predicts further migratory insertion of CO with the \mbox{--C(O)OiPr} species yielding species (4) and (8) to be also thermodynamically favorable. Subsequent energy refinement at the DFT level restores the agreement with the experimental observations. Similar reaction paths were identified for more complex catalytic model containing the molecularly defined Mn(I)-catalyst stabilized by a bidentate diamino ligand, the alkoxide base, two isopropanol solvent molecules and acetophenone as a model substrate. Simulations suggested that the nucleophilic attack of the alkoxide anion by the Mn-bound carbonyl ligand may initiate reaction paths resulting in a (partial) decoordination of the organic ligand, which can be considered the onset of catalyst deactivation~\cite{hashemi2022renegate}. These computational insights have inspired the experimental finding on the stabilizing role of Lewis acid additives mediating the nucleophilicity of the alkoxide base allowing to considerably extend the lifetime of the homogeneous Mn(I) catalysts~\cite{van2021catalysis} and the discovery of new Mn-mediated C--C coupling chemistry~\cite{van2021manganese}. The application of our reaction network analysis approach was further extended towards high-throughput computational reactivity exploration and automation identification of the classes of reactivity patterns within specific catalyst groups (Figure~\ref{fig:workflow}a)~\cite{hashemi2023hirex}. We have applied this workflow to a virtual library containing 576 Mn-pincer complexes corresponding to four distinct pincer families with varied functionalization of the ligand backbone (R1 and R2 functionalities) and Mn coordination (X) (Figure~\ref{fig:workflow}b). The 2D-MolGraphs were used to analyze the reactive trajectories, as well as to featurize and label the discovered new configurations and intermediates following the changes in the interaction patterns observed during the transformations. K-mode clustering analysis on the resulting labeled database (Figure~\ref{fig:workflow}b) has provided new insights into the reactivity of Mn(I) pincers and how it is affected by the structural modification of the ligand\unskip\break backbone. \begin{figure*} {\vspace*{-2pt}} \includegraphics{fig18} {\vspace*{-2pt}} \caption{\label{fig:workflow}(a)~A perspective workflow for the chemical space exploration and automated reactivity analysis of virtual catalyst libraries involving 2D-MolGraphs as the core component of the dynamic trajectory analysis and structure featurization. The method was used for (b)~the reactivity screening of a virtual library of 576~Mn(I) pincer complexes followed by K-mode clustering analysis of the resulting data. The figure presents top 10 clusters showing the most frequent types of interaction changes as a function of different functionalization at R1 site for each pincer class. The color coding is given at the bottom of the figure.} {\vspace*{-2pt}} \end{figure*} Specifically, the calculations indeed revealed multiple paths involving the nucleophilic attack on the carbonyl ligand and decoordination of the pincer ligand. Depending on the ligand's nature and structure, the thermodynamic favorability of such secondary conversion paths varies greatly providing thus an opportunity to tune the stability and reactivity of the transition metal catalyst, and guide the exploration of new chemical conversion\break paths. \section{Prospects and new developments in progress}\label{sec:concl} With this review, we believe that the topological 2D-MolGraphs and associated graph algorithms have been shown to be powerful tools for analyzing atomistic molecular dynamics trajectories and extracting the actual conformations sampled over time. Demonstrations were carried out on gas phase flexible molecules and inhomogeneous aqueous interfaces in the condensed phase. We have also shown the relevance of different levels of granularity to be used in the topological graphs. In particular, the coarse-grained approach of graphs of polymorphic H-bonded cycles was shown decisive in order to \mbox{rationalize} the dynamics of a flexible hexapeptide. The ``metastructures'' over which the dynamics of this peptide is built could be found only at this level of topology representation. We have also shown that the 2D-MolGraph approach can be easily coupled to global workflows that include several theoretical methods to sample conformational and reactive chemical spaces. The topological graphs were inserted in high-throughput \mbox{in silico} chemical reactive workflows in homogeneous catalysis. The conformational fingerprinting provided by the 2D-MolGraphs was decisive in several steps of these workflows. The outcome of these reactive workflows would not have been as easy and successful without the topological graphs. These works continue. Ongoing works and developments also include for instance the use of the 2D-MolGraphs to extract the knowledge of structural motifs at aqueous interfaces. For instance, it is of utmost importance to reveal the motifs formed between the sites at the surface of a solid (for instance silica oxide as in~\cite{Gaigeot_JCTC2012,Pezzotti_PCCP2019}) and the water molecules located at the interface with the solid in the Binding Interfacial Layer (BIL). These motifs are not only responsible for the spectroscopic signatures recorded at aqueous interfaces, typically by SFG (Sum Frequency Generation) spectroscopies, they are also involved in the chemical reactivity of these aqueous interfaces. Motifs and their vibrational fingerprints have already been discussed in~\cite{Pezzotti_PCCP2019} for aqueous silica interfaces. We are developing algorithms that can \mbox{automatically} recognize and classify these motifs from the 2D-MolGraphs, for instance in terms of the sizes of the H-bonded cycles formed between surface sites and water molecules and in terms of their statistical \mbox{distribution} in space within the BIL. The same algorithms will also be applied to the water molecules in the BIL of biomolecules. This is also ongoing work, using the atomistic 2D-MolGraphs and the coarse-grained polygraphs of H-bonded cycles. Naturally, databases of 2D-MolGraphs and coarse-grained polygraphs can be built up and connected to AI (Artifical Intelligence) and machine learning techniques. This is where our next step will take us. \section*{Declaration of interests} The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations. \back{} \bibliographystyle{crchim} \bibliography{crchim20240075} \refinput{crchim20240075-reference.tex} \end{document}