# Predicting redox potentials by graph-based machine learning methods

Linlin Jia, Éric Brémond, Larissa Zaida, Benoit Gaüzère, Vincent Tognetti, Laurent Joubert

► To cite this version:

Linlin Jia, Éric Brémond, Larissa Zaida, Benoit Gaüzère, Vincent Tognetti, et al. Predicting redox potentials by graph-based machine learning methods. Journal of Computational Chemistry, 2024, 45 (28), pp.2383 - 2396. 10.1002/jcc.27380. hal-05281754

HAL Id: hal-05281754

https://hal.science/hal-05281754v1

Submitted on 24 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

<figure><img src="image_3.png" /></figure>

WILEY

# Predicting redox potentials by graph-based machine learning methods

Linlin Jia<fnref n="1" /> ID | Éric Brémond<fnref n="2" /> ID | Larissa Zaida<fnref n="2" /> | Benoit Gaüzère<fnref n="3" /> ID |
Vincent Tognetti<fnref n="4" /> ID | Laurent Joubert<fnref n="4" />

**Correspondence**

Vincent Tognetti, Normandy Univ., COBRA UMR 6014 & FR 3038, Université de Rouen, INSA Rouen, CNRS, 1 rue Tesnière, 76821 Mont St Aignan Cedex, France.

Email: vincent.tognetti@univ-rouen.fr

**Abstract**

The evaluation of oxidation and reduction potentials is a pivotal task in various chemical fields. However, their accurate prediction by theoretical computations, which is a complementary task and sometimes the only alternative to experimental measurement, may be often resource-intensive and time-consuming. This paper addresses this challenge through the application of machine learning techniques, with a particular focus on graph-based methods (such as graph edit distances, graph kernels, and graph neural networks) that are reviewed to enlighten their deep links with theoretical chemistry. To this aim, we establish the ORedOx159 database, a comprehensive, homogeneous (with reference values stemming from density functional theory calculations), and reliable resource containing 318 one-electron reduction and oxidation reactions and featuring 159 large organic compounds. Subsequently, we provide an instructive overview of the good practice in machine learning and of commonly utilized machine learning models. We then assess their predictive performances on the ORedOx159 dataset through extensive analyses. Our simulations using descriptors that are computed in an almost instantaneous way result in a notable improvement in prediction accuracy, with mean absolute error (MAE) values equal to 5.6 kcal mol⁻<fnref n="1" /> for reduction and 7.2 kcal mol⁻<fnref n="1" /> for oxidation potentials, which paves a way toward efficient *in silico* design of new electrochemical systems.

**KEYWORDS**

density functional theory, graph-based machine learning methods, ORedOx159 database, Redox potential prediction

# 1 | INTRODUCTION

The experimental optimization of chemical reagents is very often a time-consuming and financially expensive task as it requires numerous tries that can also involve hazardous compounds or complex synthetic strategies. As a consequence, the exploration of the chemical space

for a given property frequently remains limited to a small number of variations, and the fine-tuning that is performed is then far from being optimal. It is thus highly desirable to have at disposal a fast screening tool that can efficiently guide the applied chemists for the selection of the best synthetic targets.

Numerical techniques are certainly suitable candidates for short-cut strategies that can be led at larger scales, provided the associated computations can treat the systems both in a reasonable time and

<footnote n="1">¹The PRG Group, Institute of Computer Science, University of Bern, Bern, Switzerland</footnote>

<footnote n="2">²Université Paris Cité, ITODYS, CNRS, Paris, France</footnote>

<footnote n="3">³LITIS, Univ Rouen Normandie, INSA Rouen Normandie, Université Le Havre Normandie, Normandie Univ, Rouen, France</footnote>

<footnote n="4">⁴Normandy Univ., COBRA UMR 6014 & FR 3038, Université de Rouen, INSA Rouen, CNRS, Mont St Aignan Cedex, France</footnote>

<footnote>A major part of this work was done when the author was at COBRA lab, France.</footnote>
with a sufficient accuracy. Advanced quantum chemistry (QC) methods, such as density functional theory (DFT) or post-Hartree-Fock methods, are indeed capable of achieving the latter point mentioned. However, their application to extensive compound searches is limited by the significant computational time required-calculations for a single molecule can span several hours to days, especially for extended systems. This limitation underscores the challenge of using these methods for broad compound exploration.

Such limitations can be alleviated by the use of machine learning (ML) methods, which usually provide predictive models that can be deployed at a large scale (then enhancing the exploration of the chemical space¹) in a rather small amount of time when they are based on features that can be evaluated in a faster way than the related full QC calculation. The importance of ML in theoretical chemistry has tre mendously increased in the last decade, and has now become ubiqui tous in the field (see for instance the reviews by Keith et al.² and Mater et al.³), generating also a lot of individual tools and software.⁴ Such booming is so exponential that it has become almost impossible to review the use of ML in chemistry in a comprehensive way, since it ranges from drug design to retro-synthesis analysis, encompassing material sciences and catalysis, or molecular dynamics and even QC itself, among other fields of application.

In this paper intended for the special issue of the Journal of Com putational Chemistry devoted to ML and artificial intelligence in theo retical and computational chemistry, we will illustrate in a didactic fashion to what extent the ML machinery can be efficiently implemen ted for the prediction of one of the most fundamental chemical prop erties, namely the redox potentials. Indeed, oxidation and reduction are elementary processes that occur in almost all chemical subfields.⁵ This is obviously the case in biochemistry since many biological pro cesses involve redox reactions (as epitomized by the cellular respira tion⁶), but also in geology (formation of blast furnaces⁷), in chemical engineering and industry (we can mention cathodic protection in gal vanized steels to fight against corrosion⁸), in new product synthesis through electrochemical catalysis,⁹ and so on.

Being able to accurately predict redox properties, both for known or new molecules, is thus of utmost importance in many areas of chemistry and also in daily life (if one thinks of batteries and solar cells for instance), and ML methods have already been applied to tackle this issue in the spirit of Quantitative Structure-Property Relation ships (QSPR) approaches. A very recent publication by Fedorov and Gryn'ova reviewed in detail this topic,¹⁰ and we thus refer the inter ested reader to this remarkable paper for an up-to-date account on ML technics for redox prediction and an expanded panorama of reported models. We will restrict ourselves to only mention here few landmark papers (the choice being of course arbitrarily and too much reduced), covering various fields (from material science to biochemis try) by Kleinová et al.,¹¹ Méndez-Hernández et al.,¹² Ghule et al.,¹³ Galuzzi et al.,¹⁴ or Bhat et al.¹⁵

Our aim is actually much less ambitious in this paper that adopts a hybrid format since it is both a short review and a research paper with new results. Hence, we will not look for a general and versatile (and even less universal) model, but, conversely, we will pick out some

specific tools—mainly belonging to graph-based methods, maybe less known in the chemical community, and we will discuss their relevant for the prediction of redox potentials. Graph theory has actually be efficiently used for the prediction of many properties belonging to dif ferent fields of chemistry, ranging from mutagenicity or toxicity in medicinal chemistry to boiling points in physical chemistry (see for instance¹⁶,¹⁷). Noteworthy, we will start from the very beginning by setting up from scratch a completely new database of organic mole cules, with DFT reference values computed by ourselves.

This departs from the common practice based on pre-existing databases. Indeed, one of the usual drawbacks of them is that they might gather reference values from various origins, sometimes with- out clear source, and they may thus suffer from a lack of homogeneity that is not without incidence on the “trustability” of the results. Con versely, we will generate here our reference values with a unique and perfectly defined, controlled, and—also an important point- reproducible computational protocol. This database will be described more in details in the next section.

Then, the computational details and the ML technics used in this study will be presented in a pedagogical way, so that it can build a bridge between the two communities (namely theoretical chemistry and data science). A particular emphasis will be put on graph-based methods since they are at the heart of our original ML methods. Then, the various ML models obtained will be presented and discussed before final conclusions.

From a methodological point of view, any ML study should adhere to the FAIR (Findable, Accessible, Interoperable, Reusable) principles.¹⁸ Artrith et al. have recently reviewed best practices in ML for chemistry¹⁹ that serve as useful guidelines for our purposes. Their checklist consists of six main points: (i) data sources should be listed, publicly available with a clear identification numbers, and with possi ble biases reported; (ii) data should be cleaned using a well-defined and discussed protocol; (iii) the methods for data representations should be clearly articulated and compared with the literature; (iv) the implementation of the ML model should be provided, and the model should be compared to baseline methods and to state-of-the-art ones; (v) a clear data split between training, validation, and testing should be implemented and clearly described; (vi) the code and workflow should be made available and allow for reproducing the reported results. Accordingly, we will illustrate all these general rules throughout this paper.

## 2 | THE OREDOX159 DATABASE

The ORedOx159 database is a collection of 318 one-electron reduc tion and oxidation reactions involving 159 large-size organic com pounds routinely used as redox indicators, or involved in the development of molecular electrochemical storage and electrochemi cal sensors.²⁰ For instance, it counts viologen derivatives which are well-known from decades to reversibly change color between violet and deep blue through reduction and oxidation,²¹ or 2,2'-bipyridiniums which are used as 'electron reservoir' for
electrochemical storage,²² or phenothiazine compounds which has been recently considered as efficient redox mediators in electrochemical sensors.²³ For each organic compound A, the database collects a one-electron reduction reaction

$$A + e^- \rightarrow A^-, \quad \Delta_r G^0_{\text{Red}} \qquad (1)$$

where $\Delta_r G^0_{\text{Red}}$ denotes the (standard) Gibbs free energy variation of the reduction process at room temperature. $\Delta_r G^0_{\text{Red}}$ values are usually computed as negative except when the oxidant A is less stable than the reducer A⁻. This last issue is related to that of negative electron affinities (in general defined for vertical processes) and can also be an artefact due to the neglecting of solvent effects (see for instance the discussion by De Proft and co-workers within a DFT context²⁴).

Complementarily, the oxidation reaction twin writes

$$A \rightarrow A^+ + e^-, \quad \Delta_r G^0_{\text{Ox}} \qquad (2)$$

with $\Delta_r G^0_{\text{Ox}}$ being the Gibbs free energy variation during the oxidation process at the same temperature.

<figure><img src="image_4.png" /><figcaption>FIGURE 1 Scheme of the main organic derivatives from which the compounds of the ORedOx159 database are derived.</figcaption></figure>

<caption>FIGURE 2 Distribution of the reduction and oxidation Gibbs free energy variations (kcal mol⁻¹) covered by the ORedOx159 database and split over the four subsets, that is, 1: phenothiazines, 2: quinones, 3: 2,2'-bipyridiniums and 4: viologens.</caption>

Let us recall that the Gibbs free energy is directly related to the standard reduction potential $E^0_{\text{Red}}$ encountered in electrochemistry. In the general case, an oxidizer is indeed reduced by z electrons to form a reduced species. These two quantities are then linked by the well-known Nernst equation

$$\Delta_r G^0_{\text{Red}} = -zF E^0_{\text{Red}}, \qquad (3)$$

where F is the Faraday constant (i.e., 96,500 C mol⁻¹). It is thus equivalent to build a predictive model on Gibbs reaction free energies or on redox potentials.

Regarding its chemical diversity, the database is composed by a large variety of compounds substituted by electron donor or acceptor groups such as aliphatic, etheroxy, keto, ester and halogen substituents, and derived from four different organic families that are clustered into the phenothiazine (1), quinone (2), 2,2'-bipyridinium (3) and viologen (4) subsets. Figure 1 provides a representation of the main core derivative of each of them, and Figure 2 depicts the distribution of their respective Gibbs free energy variations in reduction and oxidation computed in the gas phase at the PBE0/def2-TZVPP level of theory.

More precisely, the phenothiazine subset counts 28 derivatives. They are 15 carbon- and 6 aza-substituted phenothiazines as well as 4 benzothiazines and 3 sulfoxided phenothiazines. Their reduction Gibbs free energy variations span between -26.31 and 10.07 kcal mol⁻¹ while their oxidation ones vary between 140.00 and 181.55 kcal mol⁻¹.

The quinone subset composes of 37 compounds. They belong to naphto- and anthra-quinones (13 and 3 derivatives, respectively) as well as isoindole-4,7-diones (11 derivatives) and other types of quinones (10 derivatives). Their reduction (resp. oxidation) Gibbs free energy variations cover a broader energy range than phenothiazine, which spans between 82.39 (resp. 55.18) and 51.84 (resp. 213.77) kcal mol⁻¹.

The 2,2'-bipyridinium subset counts 45 derivatives. It is composed of unbridged 2,2'-bipyridiniums (3 derivatives), ethene (3 derivatives), ethane (19 derivatives), propane (6 derivatives) and butane 2,2'-bipyridiniums (3 derivatives) as well as other types of 2,2'-

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 3 Correlation plots between the Gibbs free energy variations calculated at the (left) RSX-PBE0/def2-TZVPP, (middle) PBEh-3c, and (right) B97-3c levels of theory with the PBE0/def2-TZVPP reference values.</figcaption></figure>

<table><caption>TABLE 1 Mean absolute deviations (kcal mol⁻¹) for the reduction and oxidation Gibbs free energies in the ORedOx159 database between the affine relationships of Figure 3 for RSX-PBE0/def2-TZVPP, PBEh-3c and B97-3c levels of theory, with the values obtained at the PBE0/def2-TZVPP reference level.</caption><thead><tr><th></th><th>RSX-PBE0</th><th>PBEh-3c</th><th>B97-3c</th></tr></thead><tbody><tr><td>Δ<sub>r</sub>G<sup>0</sup><sub>Red</sub></td><td>1.56</td><td>1.01</td><td>0.93</td></tr><tr><td>Δ<sub>r</sub>G<sup>0</sup><sub>Ox</sub></td><td>5.38</td><td>3.44</td><td>2.20</td></tr></tbody></table>

bipyridiniums (11 derivatives). Their reduction Gibbs free energy variations span between -228.58 and -5.67 kcal mol⁻¹, while their oxidation ones vary between 182.06 and 402.04 kcal mol⁻¹.

Finally, the viologen subset is composed of 49 derivatives, belonging to core (20 derivatives), symmetric (12 derivatives) and asymmetric (7 derivatives) substituted viologens as well as other types of viologens (10 derivatives). Their reduction (oxidation) Gibbs free energy variations cover a similar energy range as 2,2'-bipyridinium. It spans between -226.69 (185.29) and -10.59 (407.36) kcal mol⁻¹.

In order to evaluate the bias introduced by the computational protocol on the Gibbs free energies, Figure 3 compares the results obtained at the PBE0/def2-TZVPP level of theory with those obtained with other density functionals. The comparison with RSX-PBE0,²⁵,²⁶ a range-separated hybrid density functional known to provide accurate estimates for ionization potentials and electronic affinities,²⁶ depicts a quasi-ideal affine relationship. Similar conclusions can be drawn from the comparison with other composite DFT approximations such as PBEh-3c²⁷ and B97-3c,²⁸ the latter being demonstrated as an excellent cost-effective approach to predict redox potentials.²⁹ This graphical observation is also confirmed by very low mean absolute deviations (MADs) while comparing the affine predicted Gibbs free energies at B97-3c level with that of PBE0/def2-TZVPP (Table 1). The error introduced is only about 0.93 kcal mol⁻¹ in reduction and slightly larger in oxidation (2.20 kcal mol⁻¹). We thus conclude that whether the computational protocol taken as

a reference to train the ML model, its accuracy will certainly remain merely unchanged.

Please note that a more complete description of the database is provided in the Supporting Information file, and the structures and respective electronic and Gibbs free energies are accessible through the GitHub platform.³⁰

## 3 | COMPUTATIONAL DETAILS

With the collected database, we now perform the necessary computations and prediction of the redox potentials. Figure 4 exhibits the workflow of the procedure, which breaks down into five steps. Three steps during the training procedure include: (I) Computing the redox potentials for each compound in training set using DFT. The corresponding computational method is detailed at the following part of this section. (II) Constructing the descriptor representation for each compound in the training set, which is detailed in Section 4.1. (III) Training a ML model with the given descriptors and redox potentials. The list of models used in this paper is described in Section 4.2. After that, a ready-to-use model is established. When a new chemical compound with unknown redox potentials is given, we (IV) construct the descriptors for this compound the same way as in step (II), and then (V) use these descriptors as the input for the trained model, and take the output of the model as the predicted redox potential.

We first present the QC methodology to generate the reference data (I). The Gibbs free energy variations in reduction and oxidation (see Equations 1 and 2) are computed in the gas phase with density-functional theory (DFT). The structures of the 159 organic derivatives are fully optimized in their neutral ground-state as well as in their reduced and oxidized states. Their geometries are thus fully relaxed. In other words, we do not consider vertical electron removal or attachment. Each of the 477 resulting minimum energy structures is then fully characterized by a frequency computation (within the harmonic approximation) from which the thermal, vibrational and entropy contributions are computed at room temperature using standard statistical physics formulae.
<figure><img src="image_7.png" /><figcaption>FIGURE 4 The workflow for the redox potential prediction.</figcaption></figure>

To further increase the diversity of the organic systems used to train the ML models, we also introduced into the optimization process a small dataset composed by 19 various organic derivatives. Dubbed DivRedOx19, this dataset adds a collection of 38 reduction and oxidation Gibbs free energy variations to the 318 brought by ORedOx159. A more complete description of the dataset is provided in the Supporting Information file, and the structures and respective electronic and Gibbs free energies are accessible through the GitHub platform.<sup>31</sup>

The PBE0 global-hybrid density functional approximation for the exchange-correlation energy<sup>32</sup> and the def2-TZVPP basis set<sup>33</sup> are chosen as a good accuracy/cost trade-off to estimate the requested energy property for this large number of compounds. All the computations are performed with the release C.01 of the Gaussian'16 program<sup>34</sup> using a tight energy threshold criterion as well as an ultrafine integration grid.

It should be underlined that only gas phase values are here considered. Obviously, experimental data are in general related to measurements in solution. However, an accurate description of solvation effects in redox processes is far from being straightforward, and is clearly beyond the scope of this paper. The same remark applies for a detailed assessment of the chosen QC methodology. We refer the interested reader to the valuable works, among others, by the groups of Gillmore<sup>35</sup> and Grimme.<sup>29</sup>

## 4 | MACHINE LEARNING FRAMEWORK

With redox potential references generated, we now provide an in-depth elucidation of the ML framework utilized in this paper. We begin with introducing the relevant descriptors, and we subsequently delineate the ML models employed, corresponding respectively to steps (II) and (III) in Figure 4.

### 4.1 | Descriptor design

As stated in the introduction, our models will be grounded in graph theory. In a nutshell, a graph G is an ordered pair, $G = (V, E)$, of vertices (also named points or nodes) and edges (also known as links or

lines) that are unordered pairs of vertices. In a chemical context, we first construct a graph from any molecule by identifying atoms as nodes and by modeling chemical bonds by edges. It is thus obvious that the usual molecular representation of a chemical compound can be straightforwardly translated into the language of graph theory.

In order to prevent any misconceptions about molecules and graph theory, it should be however stressed that such a representation does not mean that only interactions between bonded atoms exist. Actually, a full molecular graph will be in principle a complete graph (i.e., all pairs of nodes are linked-which is not a synonym of bonded-together), and one can, in principle, evaluate the interaction energies between all nodes (for instance, using Pendás' Interacting Quantum Atoms (IQA) scheme<sup>36</sup>). However, such graphs rapidly become huge since they involve $n(n-1)/2$ edges (with n the number of atoms), so that their use would introduce a computing complexity that can be prohibitive for large databases of extended molecules. Conversely, a linear molecule (i.e., without rings or cages) exhibits (as a direct consequence of the Poincaré-Hopf relationship within the framework of Bader's Quantum Theory of Atoms-In-Molecules (QTAIM)<sup>37</sup>) only $n-1$ bonds, so that linear scaling may be possible with such graph assumption. An important point to notice is that, in principle, the node properties, if relevantly chosen, can reflect both bonded and (indirectly) non-bonded interactions.

Then, once the molecular graphs established, we consider in total three types of descriptors for graph-based models in this research. As we intend this paper to have some didactic content, we have decided to describe them in more detail with no prerequisite on chemoinformatics. They are also represented in Figure 5.

**(A) Atom + Bond Types:** This kind of descriptor considers element types (for instance, carbon, hydrogen, etc.) as node features and chemical bond types (either single, double, triple, aromatic) as edge features (see more details in type B). Since these features are discrete, they are encoded in a one-hot manner. These descriptors, while affording a simplistic description of a molecule, are widely used in various databases for the evaluation of ML models, such as TUDataset.<sup>38</sup> However, this descriptor may lack the ability to capture sufficient information for effective learning.

**(B) One-Hot Features:** Such descriptors encompass common features for both atoms and chemical bonds, each of them represented as one-hot coded node and edge features, respectively. These descriptors are inspired by the featurizers in DeepChem,<sup>39</sup> based on WeaveNet paper.<sup>40</sup> They can all be directly obtained from the RDKit<sup>41</sup> software without the need for computationally intensive QC-informed descriptors.

In practice, they are actually very often directly generated from the Simplified Molecular Input Line Entry Specification (SMILES)<sup>42</sup> molecular code (see an example at the top left corner of Figure 5) using the standard rules of chemistry based on Lewis representation, without resorting to any kind of advanced calculations. For relatively simple molecules in organic chemistry, the bond pattern predicted by the SMILES code is fully consistent with those from more elaborated approaches (for instance QTAIM already mentioned, which affords a univocal definition of chemical bonds) based on QC analysis. This is
<figure><img src="image_8.png" /><figcaption>**FIGURE 5** Descriptors used for the machine learning models.</figcaption></figure>

one of the main strengths of the SMILES code, since it can be pro-
duced and processed in an almost instantaneous manner. However, as
it does not encode any 3D information, it cannot distinguish between
the different conformers (resulting, for instance, from free rotations
around single bonds) of a given molecule.

More precisely, the node features include a 32-dimensional one-
hot vector (see middle part in Figure 5). The first 10 components
encode the chemical symbol with the following possible categories: C,
N, O, F, P, S, Cl, Br, I, and others. For instance, a carbon atom will be
coded by (1 0 0 0 0 0 0 0 0 0) while an oxygen atom will be repre-
sented by (0 0 1 0 0 0 0 0 0 0). The hybridization type is encoded by
three components depending on it is sp, sp² or sp³: a carbon atom
engaged into a double bond will hence be coded as (0 1 0).

The other atom descriptors are the electronic charge (as its evalu-
ation is not straightforward-ideally it requires a QC calculation, we
decided to disregard it by setting it to zero for all atoms), its aromatic-
ity (represented by a boolean depending whether the atom belongs to
an aromatic system or not), its so-called “Degree” (corresponding
to the coordination number), the number of hydrogen atoms linked to
the studied node (one-hot coded as a six-dimension feature), and its
possible chirality (R or S absolute configuration, once again in one-hot
encoding) in the case of an asymmetric carbon atom. This representa-
tion is thus able to address stereochemistry issues.

The edge features consist of an 11-dimensional one-hot vector
(see bottom part of Figure 5) representing bond type (four categories:
single, double, triple and aromatic), whether the atoms at the end of
the bond share the same ring (no matter the type) and whether the

bond belongs to a conjugated system (both encoded as a boolean),
and the possible stereoisomerism of the bond in a five-dimension
feature.

Notice that some structure information is implicitly embedded
into the node features through “Aromatic”, “Degree”, and “Chirality”.
This may benefit the models that can not deal with edge features or
rings. We emphasize that these features are easy-to-get and non-
expensive, thus exhibiting how much ML models can benefit from a
short amount of effort of descriptor engineering.

(C) **One-Hot Features on Nodes + Distances on Bonds:** This descriptor employs the same node features as the descriptor (B), but replaces edge features with the Euclidean distances (here expressed in Angström) between the nuclei of the corresponding bonded atoms (see bottom right corner in Figure 5). We utilize distances to maintain a 3D invariant system (that-is-to-say the values remain the same when the molecule is translated and rotated as a whole), making it applicable to any models with the ability to handle continuous edge features (at variance with the previous bond descriptions that used only discrete descriptors). It can be remarked that internuclear distances are also the main ingredient of the celebrated Coulomb matrices (CMs)⁴³ that enter many ML models in chemistry. However, it should be noticed that distances, here (at variance with that is currently done when using CMs), are not computed for non-bonded pair of atoms even if they can add additional information on the geometry adopted by the molecule.

To obtain these distances, we search for the most stable con-
former for all compounds using the freely available Balloon⁴⁴ software
<figure><img src="image_9.png" /><figcaption>FIGURE 6 Machine learning models used for the prediction.</figcaption></figure>

with the MMFF94 molecular mechanics force field, starting from the SMILES code, with the following options: -nconfs 20 -nGenerations 300 -rebuildGeometry. We then compute the distance from the corresponding Cartesian coordinates of the nuclei. This computation, though probably less accurate than a QC geometry optimization, has however an almost negligible computational cost, fulfilling one of the targets of our ML modelling.

For the same reason, no QC descriptors are used at all. Obviously, some of them would be certainly useful for an accurate prediction of redox properties. For instance, the energy of the Highest Occupied Molecular Orbital (HOMO) energy is known, in Kohn-Sham DFT, to be exactly equal (this is due to the asymptotic form of the electron density in the exponential tail) to the opposite of the electronic component (thus without thermodynamic contributions) of the vertical (i.e., without geometry relaxation) ionization energy if the exact exchange-correlation functional is used (which is unfortunately not known in analytical form). It is also known that the energy of the lowest unoccupied molecular orbital (LUMO), while differing in principle from the opposite of the vertical electron affinity (as a consequence of the discontinuity of the exchange-correlation potential), is nevertheless linked to electron capture (see for instance a detailed discussion by Baerends⁴⁵).

These descriptors, which have already been used in the literature to build predictive models, unfortunately require the computationally demanding task of a DFT calculation. This is exactly what we would like to avoid thanks to a well tailored ML model based on descriptors that can be computed at an almost zero cost.

## 4.2 | Machine learning models

This section outlines the machine learning models utilized in this study, as depicted in Figure 6. Initially, we employ seven vector-based machine learning models as baselines. Subsequently, we introduce three categories of graph-based models, namely, five Graph Kernels (GKs), one Graph Edit Distance (GED) model, and six Graph Neural Networks (GNNs), making a total of twelve models. These models are systematically evaluated on their performances for the graph prediction challenges.

**Vector-based Models:** First, we select seven models that process vectors inputs as baselines.⁴⁶ These models include Linear Regression (LR), Gaussian Kernel Ridge Regression (GKRR), Support Vector Regression (SVR), Gaussian Process Regression (GPR), Random Forest (RF), Gradient-boosted Decision Trees (GBDT), and K-Nearest Neighbor Regression (KNN). Given the vector-focused nature of these models, a graph-level pooling step is required beforehand to derive fixed size vectors from varying size graphs. For this purpose, we adopted the strategy outlined in Reference 47, wherein statistics are computed for each feature over all nodes within each graph (stats). Indeed, the number of nodes varies along the molecular dataset, so that the vector gathering the node features would have, in general a different size for one chemical system to another one, precluding the use of standard ML methods. To avoid this, for a given feature, we compute the minimal, maximal, mean and standard deviation values on all nodes, resulting in four statistical descriptors that can be evaluated for all molecules.
Additionally, we introduced a novel pooling approach for the vector-based models, where we quantified the frequency of each feature. As discussed in Section 4.1, the node features, which encompass Atom and Bond Types (A) and One-Hot Features (B), are encoded using a one-hot representation in this paper. When two feature vectors are summed, it results in a vector where each element corresponds to the total count of that feature. To illustrate this concept, consider a simple two-node toy graph with feature vectors (1 0 0 1 0) and (0 1 0 1 0) assigned to the two nodes. After applying our pooling approach, the graph-level representation becomes (1 1 0 2 0), obtained through element-wise summation of the two node feature vectors.

This means that features at the first, second, and fourth positions occur once, once, and twice, respectively, in the entire graph, while the other features do not appear. For a visual representation, refer to Figure 6 (I). Our ML experiments show that this pooling strategy outperforms the one proposed in Reference 47 in most cases, as illustrated in Section 5. It is important to note that these models primarily emphasize node features and that, during the experiments, edge features were disregarded.

**Graph Kernels:** Graph kernels compute a similarity measure between graphs by implicitly mapping them into high-dimensional spaces where the inner product is computed. Such a measure can be defined as a symmetric, positive semi-definite function $k: G \times G \to \mathbb{R}$ on a space of graphs $G$, where a map $\phi: G \to \mathcal{H}$ exists into a Hilbert space $\mathcal{H}$, so that

$$k(G_1, G_2) = \langle \phi(G_1), \phi(G_2) \rangle_{\mathcal{H}} \quad (4)$$

for all pairs of graphs $G_1, G_2 \in \mathcal{G}$. Here $\langle \cdot, \cdot \rangle_{\mathcal{H}}$ is the inner product in $\mathcal{H}$, as shown in Figure 6 (II). The results are then processed by a kernel machine for final prediction, specifically Kernel Ridge Regression (KRR)$^{46}$ in this paper. By altering the design of the mapping, various graph kernels can be constructed. A commonly design strategy is inherited from R-convolution kernels, which measures similarity between two objects, by measuring the similarities between their substructures.$^{48}$

Normally, a series of sub-kernels are first constructed between pairs of sub-structures of two graphs, then a graph-level kernel is established upon the summation of these sub-kernels. In our paper, we examine a set of baseline graph kernels, including Shortest Path Kernel (SP),$^{49}$ Structural Shortest Path Kernel (SSP),$^{50}$ Path Kernel,$^{16}$ Treelet Kernel,$^{51}$ and Weisfeiler-Lehman Subtree Kernel (WLSubtree),$^{52}$ each of them named after the sub-structures from which the kernel is built. These kernels are able to tackle linear and non-linear sub-structures. Detailed descriptions and comparisons of graph kernels can be found in References 53–56.

Graph kernels are able to work directly with the graph structures, while tackling the similarity measure on high or infinite dimension spaces, thanks to the kernel trick.$^{57}$ This allows bypassing acquiring fixed embeddings, extends their expressiveness, and flexibility of design and integration of prior-knowledge. The proper design can lead to the invariance to graph isomorphism as well. The other size of the

coin minted their lack of the ability of automatic learning of the representation, especially compared to graph neural networks. The choice of the appropriate kernel is thus crucial for achieving good performance.

**Graph Edit Distances:** GEDs measure the dissimilarity between two graphs by evaluating the amount of distortion required to transform one graph into another. Various basic edit operations, each assigned with a cost, are used to measure this distortion, and the minimum total cost represents the GED between the two graphs. Commonly used edit operations include an insertion, removal, or substitution of vertices or edges. Figure 6 (III) shows a simple instance of an edit procedure. In this basic example, the first operation is the removal of the bottom vertex, followed by the removal of the bottom edge.

Then an edit path $\pi$, which in Figure 6 goes from graph $G$ to graph $G'$, can be defined as the sequence of these operations, whose total cost $c$ is the sum of the costs of all edit operations in the path, namely $c(\pi) = \sum_{e \in \pi} c(e)$, where $e$ denotes an elementary operation. Then, the GED between two graphs $G_1$ and $G_2$ is defined as the minimum cost associated with any possible edit path:

$$ged(G_1, G_2) = \min_{\pi \in \Pi(G_1, G_2)} c(\pi), \quad (5)$$

where $\Pi(G_1, G_2)$ denotes the set of all possible edit paths from $G_1$ to $G_2$. We mention here that we only consider in this paper a constant value for the cost (i.e., it is independent of the labels of nodes and edges) of every of the six elementary operations.

Since computing GED is an NP-hard problem,$^{58}$ approximation methods have been proposed to estimate it. In this paper, we use a milestone method, called Bipartite,$^{59}$ where the edit costs are optimized using a method from Reference 60. This approach boils down the quadratic assignment problem (QAP) to a linear one by finding a suboptimal edit path only based on local structural information. Edit costs are optimized in a bi-alternate optimization scheme. This approach constitutes a good trade off between computational time and accuracy. More accurate methods based on the QAP have been proposed,$^{61,62}$ at the cost of an higher computational time. We recommend in-depth papers$^{63,64}$ on GEDs to interested readers.

As GEDs are a distance measure, we combine them with the k-Nearest Neighbor regressor (KNN)$^{65}$ for the final prediction. To predict the output value of an example, this approach simply consists in averaging the target values associated to its k nearest neighbors according to the GED, with the hypothesis that the target values of similar examples are close to each other. The contribution of each neighbor can be weighted according to distance values.

The advantages of GEDs inherit from their flexible and explicit application of edit operations, which allows for capturing complex structural differences between graphs, as well as an explicit demonstration of the modification process of graphs. This latter merit helps establish the explainability of the model, which is conductive to their application on the generative tasks. However, the semantic meaning of graph elements may not be captured by GEDs. Moreover, as a NP-hard problem, calculating GEDs can be computationally expensive,
especially for large or dense graphs. The lack of providing a continuous space to interpolate between graphs and invariance to graph size may also limit their applications to certain tasks.

**Graph Neural Networks (GNNs):** GNNs are a class of neural networks that specialize in learning from graph-structured data. They often utilize a message-passing scheme, wherein information from nodes, edges, and the overall graph structure, are aggregated from the neighbors of each node through a series of graph convolutional operations, as shown in Figure 6 (IV). With denoting $\mathbf{x}_i^{(k-1)} \in \mathbb{R}^F$ node features of node i in layer $(k-1)$, and $\mathbf{e}_{j,i} \in \mathbb{R}^D$ denoting (optional) edge features from node j to node i, the updating of node features using message passing graph neural networks can be described as

$$ \mathbf{x}_j^{(k)} = \gamma^{(k)}(\mathbf{x}_j^{(k-1)}, \bigoplus_{j \in \mathcal{N}(i)} \phi^{(k)}(\mathbf{x}_i^{(k-1)}, \mathbf{x}_j^{(k-1)}, \mathbf{e}_{j,i})) \quad (6) $$

where $\oplus$ denotes a differentiable, permutation invariant function (e. g., sum, mean or max), and $\phi$ denote differentiable functions such as Multi Layer Perceptrons (MLPs). $\gamma^{(k)}$ can be thought of as the node's update or processing function at that particular layer, which can be as well implemented as a neural network layer. This formula captures the message-passing operation of GNNs. At each layer, information is propagated from neighboring nodes to the central node i.⁶⁶ Following this message-passing process, a pooling operation (such as sum, mean, or max) over node features is often applied to obtain a graph-level representation; a MLP can then be used for graph-level predictions. We refer to References 67–70 readers interested in GNNs.

Several representative GNN architectures have been considered for study. The Message Passing Neural Network (MPNN)⁷¹ employs message passing and aggregation mechanisms in combination with an edge network to capture relational information. In contrast, the Graph Convolutional Network (GCN)⁷² utilizes convolution-like operations with weight sharing for neighbor node information aggregation. The Deep Graph Convolutional Neural Network (DGCNN)⁷³ distinguishes itself through the use of edge features and sort pooling to learn graph embeddings invariant to permutation of the nodes. The Graph Isomorphism Network (GIN)⁷⁴ combines a graph isomorphism test with learnable functions to integrate global and local information. On the other hand, the Graph Attention Network (GAT)⁷⁵ employs attention mechanisms to selectively weight neighboring node contributions.

Finally, the Unified Message Passing Framework (UniMP)⁷⁶ integrates message passing and transformer-based architectures, harnessing both relational and self-attention mechanisms for comprehensive graph representation learning. These diverse models offer distinct approaches to address our problem.

One significant advantage of GNN models lies in their time complexity during inference, setting them apart from GEDs and graph kernels. In the case of GEDs and graph kernels, it is typically necessary to perform metric or Gram matrix computations, involving the calculation of similarity values between each data point in the training set and the data point to be inferred. This matrix serves as the inner product in kernel machines and aids in finding the K nearest neighbors in K-Nearest Neighbors algorithms. The catch is that the computational

time associated with this matrix scales linearly with the size of the training set, rendering it prohibitively time-consuming for large datasets. As a result, the scalability of these models becomes constrained.

In contrast, GNN models usually feature a linear or nearly linear time complexity during inference. This efficiency is intricately linked to the model's design, particularly the scale of weight parameters. In practice, GNN computations predominantly involve basic matrix operations like multiplication and addition. Consequently, GNN models excel in efficiently handling graphs of varying sizes, making them exceptionally well-suited for real-world scenarios. These scenarios often involve graph data that ranges from small and simple to large and intricate, then underlining the substantial benefits of GNNs in real-world applications.

More advantages are offered by GNNs, notably in their ability to automatically learn representations, which alleviates the need for explicit feature engineering. By aggregating and propagating attributes associated with nodes and edges within the graph, the incorporation of both local and global information is facilitated for a comprehensive understanding of the graph's context. Additionally, GNNs allow for knowledge transfer through a pre-training and fine-tuning strategy. However, they come with certain limitations. GNNs often require substantial amounts of training data and memory resources, which can be restrictive in real-world applications.

The lack of proper interpretability and explainability poses challenges, potentially diminishing trust among experts. Furthermore, GNNs inherit unresolved fundamental issues, such as over-smoothing, where nodes tend to lose their distinctiveness after a certain number of GNN layers, limiting the feasibility of designing deep GNN architectures. These complexities need further investigation and innovation to fully harness GNNs' potential in practical settings.

## 5 | RESULTS AND DISCUSSION

We first mix the ORedOx159 and DivRedOx19 datasets, dubbed Mix-RedOx178 in the rest of the paper, to increase the diversity of the considered organic systems. Then performance evaluations are carried out on ten different random splits. Each split is partitioned into 80%, 10%, and 10% for training, validation, and testing, respectively. The hyper-parameters for each model are tuned on the training sets separately for every split using a grid-search approach. Consequently, 10 unique instances, each with a different set of hyper-parameters, are selected for each model, corresponding to one instance per split. We calculate the Average Mean Absolute Errors (MAEs) over the 10 splits as the final results, estimated on the respective test sets.

Table 2 shows the computing resource settings for the experiments, in which GCCore⁷⁷ and Cython⁷⁸ are for GED computations, and Pytorch⁷⁹ and Pytorch Geometric⁶⁶ packages are used to build GNN models. A more detailed description of the experimental settings, including the tuning strategy, tuning ranges, and explanation of the hyper-parameters, and the scales of weight parameters of the GNN models, can be found in the Supporting Information file. For the sake of reproducibility, the code is available through the GitHub platform.⁸⁰
Table 3 presents the MAE values obtained for each ML experiment. The ± sign indicates the 95% confidence interval computed over the 10 repetitions. Unsupported descriptors for each model are denoted by “-”. Notably, the best descriptor for each model is highlighted in bold, the optimum result for each descriptor across all models is marked in green, and the superior results across all experiments are underlined. In the case of vector-based models, we focus on the results obtained through count pooling as it consistently outperforms statistical pooling. This count pooling thus represents an interesting improvement over the methods previously used

by us to deal with the non-constant number of nodes along the whole dataset. An exception is observed for the GPR model with the One-Hot descriptor, where statistical pooling yields better results, and thus we present both pooling methods in this case.

When comparing different descriptors, a clear trend emerges with the One-Hot descriptor prevailing in most models: 16 out of 20 for both reduction and oxidation Gibbs free energies. In the context of graph-based models, this dominance continues as 8 out of 12 models for the reduction targets and 9 out of 12 for the oxidation targets. Importantly, the performance gap among descriptors is significant on a model-wise basis, which often surpasses the differences in performance between models. This phenomenon is particularly pronounced for vector-based models, likely due to their limitation in handling only node features. The One-Hot descriptor, encoding implicit structural information, proves highly beneficial in this context.

A similar pattern is observed for most GNN models. Conversely, the Distance descriptor is evaluated only for models supporting continuous edge features. Its performance closely aligns with that of the One-Hot descriptor, sometimes surpassing the latter. These observations underscore the critical role of descriptor engineering in performance.

Further insights emerge when comparing different models. Notably, all optimal results are obtained by graph-based models,

<table><caption>TABLE 2 Computing resource settings for experiments.</caption><thead><tr><th>Environments</th><th>Settings</th></tr></thead><tbody><tr><td>CPU</td><td>128x AMD EPYC 7742<br>64-Core Processor</td></tr><tr><td>Operating system</td><td>CentOS 7.9.2009.x86_64</td></tr><tr><td>Python version</td><td>3.8.6</td></tr><tr><td>GCCcore version</td><td>12.2.0</td></tr><tr><td>Cython version</td><td>0.29.21</td></tr><tr><td>PyTorch version</td><td>2.0.1+cu118</td></tr><tr><td>PyTorch Geometric version</td><td>2.3.1</td></tr></tbody></table>

<table><caption>TABLE 3 The prediction MAE (in kcal/mol).</caption><thead><tr><th rowspan="2">Categories</th><th rowspan="2">Models</th><th colspan="3">Δ<sub>r</sub>G<sup>0</sup><sub>Red</sub></th><th colspan="3">Δ<sub>r</sub>G<sup>0</sup><sub>Ox</sub></th></tr><tr><th>AB types</th><th>One-Hot</th><th>Dis</th><th>AB types</th><th>One-Hot</th><th>Dis</th></tr></thead><tbody><tr><td rowspan="7">Vector-Based</td><td>LR</td><td>56.0±10.8</td><td><strong>29.4±2.7</strong></td><td>—</td><td>56.1±11.1</td><td><strong>31.0±2.5</strong></td><td>—</td></tr><tr><td>GKRR</td><td>28.9±7.5</td><td><strong>14.4±5.0</strong></td><td>—</td><td>30.5±10.4</td><td><strong>17.5±4.4</strong></td><td>—</td></tr><tr><td>SVR</td><td>26.0±9.2</td><td><strong>14.1±3.7</strong></td><td>—</td><td>31.4±10.3</td><td><strong>17.2±5.2</strong></td><td>—</td></tr><tr><td>GPR (stats)</td><td>42.3±6.5</td><td><strong>21.0±5.5</strong></td><td>—</td><td>43.2±8.8</td><td>22.0±3.0</td><td>—</td></tr><tr><td>GPR (count)</td><td>36.8±7.3</td><td><strong>26.2±4.6</strong></td><td>—</td><td>37.6±8.8</td><td><strong>29.6±5.1</strong></td><td>—</td></tr><tr><td>RF</td><td>19.7±5.0</td><td><strong>14.3±4.6</strong></td><td>—</td><td>22.3±5.4</td><td><strong>15.2±4.0</strong></td><td>—</td></tr><tr><td>GBDT</td><td>19.6±5.1</td><td><strong>11.9±4.1</strong></td><td>—</td><td>20.4±4.9</td><td><strong>13.7±3.2</strong></td><td>—</td></tr><tr><td>KNN</td><td></td><td>25.3±9.9</td><td><strong>10.2±3.9</strong></td><td>—</td><td>26.2±8.2</td><td><strong>11.3±5.4</strong></td><td>—</td></tr><tr><td rowspan="5">Graph Kernels</td><td>SP</td><td>29.6±14.4</td><td><strong>12.5±2.1</strong></td><td>—</td><td>31.9±13.9</td><td><strong>14.0±2.1</strong></td><td>—</td></tr><tr><td>SSP</td><td><strong>9.5±2.3</strong></td><td>14.7±2.8</td><td>14.4±3.2</td><td><strong>11.9±2.7</strong></td><td>16.6±2.8</td><td>16.9±3.4</td></tr><tr><td>Path</td><td><strong>9.7±2.1</strong></td><td>10.7±2.4</td><td>—</td><td><strong>12.8±2.7</strong></td><td><strong>12.8±2.3</strong></td><td>—</td></tr><tr><td>Treelet</td><td><strong>10.3±2.7</strong></td><td>13.8±3.2</td><td>—</td><td><strong>13.4±3.1</strong></td><td>17.3±3.5</td><td>—</td></tr><tr><td>WLSubtree</td><td>12.4±2.1</td><td><strong>12.0±1.8</strong></td><td>—</td><td><strong>12.0±1.8</strong></td><td>13.6±2.2</td><td>—</td></tr><tr><td>GEDs</td><td>Bipartite (fitted)</td><td>11.9±4.2</td><td><strong>8.0±1.8</strong></td><td>9.5±4.0</td><td>14.3±5.1</td><td><strong>12.1±2.7</strong></td><td>13.8±4.9</td></tr><tr><td rowspan="6">GNNs</td><td>MPNN</td><td>8.2±3.7</td><td><strong>5.8±1.9</strong></td><td><u>5.6±0.6</u></td><td><strong>8.7±3.2</strong></td><td><u>7.2±2.5</u></td><td>10.7±3.3</td></tr><tr><td>GCN</td><td>21.5±7.6</td><td><strong>8.4±3.2</strong></td><td>—</td><td>24.5±8.4</td><td><strong>8.8±2.0</strong></td><td>—</td></tr><tr><td>DGCNN</td><td>17.7±6.2</td><td><strong>8.8±3.2</strong></td><td>—</td><td>22.1±6.2</td><td><strong>13.5±4.8</strong></td><td>—</td></tr><tr><td>GIN</td><td>26.9±9.0</td><td><strong>13.7±5.0</strong></td><td>—</td><td>24.7±9.7</td><td><strong>13.0±2.9</strong></td><td>—</td></tr><tr><td>GAT</td><td>11.5±4.9</td><td><strong>6.9±1.4</strong></td><td>8.0±2.3</td><td>11.1±6.6</td><td><strong>8.5±1.5</strong></td><td>9.9±2.2</td></tr><tr><td>UniMP</td><td><strong>7.3±1.5</strong></td><td><strong>7.4±3.0</strong></td><td>7.9±1.0</td><td>18.0±18.3</td><td><strong>9.1±2.6</strong></td><td>20.3±24.5</td></tr></tbody></table>

Note: The best descriptor for each model is highlighted in bold, the optimum result for each descriptor across all models is marked in green, and the superior results across all experiments are underlined.
particularly GNNs when using the One-Hot and the Dis descriptors. The best performances across all experiments for both targets are achieved by MPNN, on the Dis descriptor for the reduction and on the One-Hot for the oxidation, underlining the flexible capabilities of graph-based models and the representational power of GNNs. Meanwhile, the GAT model performs quite well on all three descriptors. This suggests the intriguing potential of attention-based models in capturing core information for graphs with limited features, even on smaller datasets. These insights emphasize the importance of finding the correct synergy between descriptors and models, indicating that their successful integration is essential for optimal performance.

On the whole, with the chosen low-level and computational-friendly descriptors, the best ML model achieves a prediction MAE of 5.6 kcal mol⁻¹ for reduction and 7.2 kcal mol⁻¹ for oxidation on the MixRedOx178 database. It corresponds to a MAE ranging between 0.2 and 0.3 V on the potential, an error which is in line with the state-of-the-art approach recently developed in computational chemistry.²⁹

<figure><img src="image_10.png" /><figcaption>FIGURE 7 Potentials computed by DFT versus the ones predicted via the best descriptor and model.</figcaption></figure>

As previously described, our experimental setup conducts performance evaluations using 10 random splits of the dataset, with each split comprising 178 compounds and designating 10% of these for testing. This arrangement results in 18 test data-points for each split, leading to an aggregate of 180 test data-points from the MixRedOx178 dataset across all the splits. Subsequently, we generated plots for these test data, comparing the Gibbs free energies calculated by DFT with their predicted counterparts, as shown in Figure 7. These comparisons were made using the best-performing prediction system, which includes the MPNN model, as well as the Dis descriptor and the One-Hot descriptor for reduction and oxidation tasks, respectively. The figure displays high correlations between the DFT-calculated and the predicted values. Specifically, the R² scores amount to 0.992 for reduction and 0.988 for oxidation. These substantial correlation coefficients validate the accuracy of our predictive model for these tasks and affirm its predictive abilities. This level of agreement reflects the model's suitability for practical applications in predicting reduction and oxidation potentials.

Lastly, we briefly discuss the computational effort associated to these ML procedures. Figure 8 illustrates the reference time required for each datapoint or compound. The values are presented in seconds and log-scaled with a base of 10. In this representation, faster reference times, characterized by smaller values, are indicated by the blue color, while slower times are depicted in red. Among the evaluated models, GNN models demonstrate superior efficiency in terms of reference times compared to graph kernel and GED models. For GNN models, reference times can be as low as 10⁻⁴ s, as seen in the case of GCN, DGCNN, and GIN. The most accurate model, MPNN, achieves reference times of approximately 0.001 s for reduction and 0.01 s for oxidation. These values represent, as expected, a dramatic

<figure><img src="image_11.png" /><figcaption>FIGURE 8 Reference time per datapoint for graph-based models, in log-scaled with base 10, and in seconds.</figcaption></figure>
improvement in efficiency compared to quantum chemistry-based methods, for which minutes or hours are needed for a molecule.

## 6 | CONCLUSION AND FUTURE WORK

In this study, we proposed a comprehensive, homogeneous, and reliable database ORedOx159 (based on density functional theory calculations), which collects one-electron reduction and oxidation reactions involving large-size organic compounds of experimental interest. We then applied a series of machine learning models, with a particular emphasis on graph theory-based methods, on this dataset, extended by 19 various organic derivatives, to predict the reduction and oxidation Gibbs free energies using several descriptors computed from the structures of the compounds obtained at the cheap molecular mechanics level.

The analyses of the results suggest that the use of (almost) instantly computed descriptors has achieved notably high prediction accuracy, incorporated with the graph-based models. The Message Passing Neural Network (MPNN) among all models achieves the best results: mean absolute errors of 5.6 kcal/mol for reduction and 7.2 kcal/mol for oxidation potentials. Particularly, the One-Hot descriptor consistently enhances model performance, prevailing all other descriptors. These outcomes not only demonstrate the crucial impact of thoughtful descriptor design, but also illustrate the effectiveness of tailored model-descriptor combinations.

Looking ahead, future research work will include several dimensions. Firstly, from the chemical point of view, the next important step will be to (i) include to our model solvation effects, and (ii) extend the training database to other chemical families of high interest in electrochemistry. Then, in comparison to tabulated reference potentials, our research output will become a robust built-in package to fast and accurately predict redox potentials of organic compounds in solutions, as it is for instance the case for nuclear magnetic resonance spectrum simulations.

Secondly, state-of-the-art machine learning models may serve as better prediction tools, especially the ones tailored for our problem. For example, Graph Neural Networks and transformer-like models, particularly focusing on leveraging 3D coordinates may unlock new capabilities and insights. Other machine learning strategies, such as pre-training, may also unlock the potential of better prediction. Thirdly, we aim at exploring more suitable descriptors that can better balance the representation ability and the acquiring time complexity. Lastly, we plan to broaden our evaluation criteria to incorporate additional important metrics such as robustness, interpretability and explainability, and overall model trustworthiness.

Finally, from a more general perspective, we are convinced that significant progress can be made mainly if the two involved scientific communities (i.e., theoretical chemistry researchers and data scientists) enter a fruitful dialogue, and manage not only to share their tools but also to build a common language.

## ACKNOWLEDGMENTS

E.B. and L.Z. thank Dr. P.P. Lainé (ITODYS lab, CNRS) for fruitful discussions regarding the construction of the ORedOx159 database.

E.B. gratefully acknowledges ANR (Agence Nationale de la Recherche) for the financial support of this work through the MoMoPlasm project (ANR-21-CE29-0003). He thanks also ANR and CGI (Commissariat à l'Investissement d'Avenir) for their financial support of this research through Labex SEAM (Science and Engineering for Advanced Materials and devices) ANR-10-LABX-096, ANR-18-IDEX-0001, and acknowledges TGGC (Très Grand Centre de Calcul du CEA) for computational resources allocation (Grant No. AD010810359R1).

This work has also been partially supported by University of Rouen Normandy, INSA Rouen Normandy, the "Centre National de la Recherche Scientifique" (CNRS), the European Regional Development Fund (ERDF), Labex SynOrg (ANR-11-LABX-0029), Carnot Institut I2C, the graduate school for research XL-Chem (ANR-18-EURE-0020 XL CHEM), and the "Région Normandie". The Centre Régional Informatique et d'Applications Numériques de Normandie (CRIANN) is acknowledged for providing access to computational resources to L.J., V.T. and L.J.

## DATA AVAILABILITY STATEMENT

The data that supports the findings of this study are available in the supplementary material of this article

## ORCID

Linlin Jia https://orcid.org/0000-0002-3834-1498

Éric Brémond https://orcid.org/0000-0002-8646-9365

Benoit Gaüzère https://orcid.org/0000-0001-9980-2641

Vincent Tognetti https://orcid.org/0000-0002-0292-3509

## REFERENCES

[1] C. W. Coley, *Trends Chem. 2021*, 3, 133.

[2] J. A. Keith, V. Vassilev-Galindo, B. Cheng, S. Chmiela, M. Gastegger, K.-R. Müller, A. Tkatchenko, *Chem. Rev.* **2021**, 121, 9816.

[3] A. C. Mater, M. L. Coote, *J. Chem. Inf. Model.* **2019**, 59, 2545.

[4] A. Hagg, K. N. Kirschner, *J. Chem. Inf. Model.* **2023**, 63, 4505.

[5] J. Schüring, H. D. Schulz, W. R. Fischer, J. Böttcher, W. H. Duijnisveld, *Redox: Fundamentals, Processes and Applications*, Springer-Verlag, Berlin and Heidelberg GmbH **1999**.

[6] L. Stryer, *Biochemistry*, (fourth ed.) ed., Freeman and Company, New York **1995**.

[7] F. Oeters, M. Ottow, H. Meiler, H. B. Lüngen, M. Koltermann, A. Buhr, J.-i. Yagi, L. Formanek, F. Rose, J. Flickenschild, R. Hauk, R. Steffen, R. Skroch, G. Mayer-Schwinning, H.-L. Bünnagel, H.-G. Hoff, *Ullmann's Encyclopedia of Industrial Chemistry*, John Wiley & Sons, Ltd, Hoboken **2006**.

[8] W. Baeckmann, W. Schwenck, W. Prinz, *Handbook of Cathodic Corrosion Protection*, 3rd ed., Elsevier, Houston **1997**.

[9] Y. Holade, K. Servat, S. Tingry, T. W. Napporn, H. Remita, D. Cornu, K. B. Kokoh, *ChemPhysChem* **2017**, 18, 2573.

[10] R. Fedorov, G. Gryn'ova, *J. Chem. Theory Comput.* **2023**, 19, 4796.

[11] M. Kleinova, M. Hewitt, V. Brezova, J. C. Madden, M. T. D. Cronin, *Gen. Phys. Biophys.* **2006**, 26, 97.

[12] D. D. Méndez-Hernández, J. G. Gillmore, L. A. Montano, D. Gust, T. A. Moore, A. L. Moore, V. Mujica, *J. Phys. Org. Chem.* **2015**, 28, 320.
[13] S. Ghule, S. R. Dash, S. Bagchi, K. Joshi, K. Vanka, ACS Omega **2022**, *7*, 11742.

[14] B. G. Galuzzi, A. Mirarchi, E. L. Viganò, L. De Gioia, C. Damiani, F. Arrigoni, *J. Chem. Inf. Model.* **2022**, *62*, 4748.

[15] V. Bhat, P. Sornberger, B. S. S. Pokuri, R. Duke, B. Ganapathysubramanian, C. Risko, *Chem. Sci.* **2023**, *14*, 203.

[16] L. Ralaivola, S. J. Swamidass, H. Saigo, P. Baldi, *Neural Netw* **2005**, *18*, 1093.

[17] B. Gaüzère, L. Brun, D. Villemin, *Pattern Recognit Lett* **2012**, *33*, 2038.

[18] M. D. Wilkinson, M. Dumontier, I. J. Aalbersberg, G. Appleton, M. Axton, A. Baak, N. Blomberg, J.-W. Boiten, L. B. da Silva Santos, P. E. Bourne, J. Bouwman, A. J. Brookes, T. Clark, M. Crosas, I. Dillo, O. Dumon, S. Edmunds, C. T. Evelo, R. Finkers, A. Gonzalez-Beltran, A. J. Gray, P. Groth, C. Goble, J. S. Grethe, J. Heringa, P. A. Hoen, R. Hooft, T. Kuhn, R. Kok, J. Kok, S. J. Lusher, M. E. Martone, A. Mons, A. L. Packer, B. Persson, P. Rocca-Serra, M. Roos, R. van Schaik, S.-A. Sansone, E. Schultes, T. Sengstag, T. Slater, G. Strawn, M. A. Swertz, M. Thompson, J. van der Lei, E. van Mulligen, J. Velterop, A. Waagmeester, P. Wittenburg, K. Wolstencroft, J. Zhao, B. Mons, *Sci Data* **2016**, *3*, 160018.

[19] N. Artrith, K. T. Butler, F.-X. Coudert, S. Han, O. Isayev, A. Jain, A. Walsh, *Nat. Chem.* **2021**, *13*, 505.

[20] P. Wardman, *J. Phys. Chem. Ref. Data* **1989**, *18*, 1637.

[21] L. Michaelis, E. S. Hill, *J. Gen. Phys.* **1933**, *16*, 859.

[22] M. Hromadová, P. P. Lainé, *Curr. Opin. Electrochem.* **2022**, *34*, 100996.

[23] M. M. Barsan, M. E. Ghica, C. M. Brett, *Anal. Chim. Acta* **2015**, *881*, 1.

[24] F. De Profit, N. Sablon, D. J. Tozer, P. Geerlings, *Faraday Discuss.* **2007**, *135*, 151.

[25] E. Brémond, A. J. Pérez-Jiménez, J. C. Sancho-García, C. Adamo, *J. Chem. Phys.* **2019**, *150*, 201102.

[26] E. Brémond, A. J. Pérez-Jiménez, J. C. Sancho-García, C. Adamo, *J. Chem. Phys.* **2020**, *152*, 244124.

[27] S. Grimme, J. G. Brandenburg, C. Bannwarth, A. Hansen, *J. Chem. Phys.* **2015**, *143*, 054107.

[28] J. G. Brandenburg, C. Bannwarth, A. Hansen, S. Grimme, *J. Chem. Phys.* **2018**, *148*, 064104.

[29] H. Neugebauer, F. Bohle, M. Bursch, A. Hansen, S. Grimme, *J. Phys. Chem. A* **2020**, *124*, 7166.

[30] Github, Structures and energies the oredox159 database are freely accessible following the link. https://github.com/ANRMoMoPlasm/ ORedOx159

[31] Github, Structures and energies the divredox19 database are freely accessible following the link. https://github.com/ANRMoMoPlasm/ DivRedOx19

[32] C. Adamo, V. Barone, *J. Chem. Phys.* **1999**, *110*, 6158.

[33] F. Weigend, R. Ahlrichs, *Phys. Chem. Chem. Phys.* **2005**, *7*, 3297.

[34] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, D. J. Fox, *Gaussian 6 Revision C.01*, Gaussian Inc, Wallingford CT **2016**.

[35] E. J. Lynch, A. L. Speelman, B. A. Curry, C. S. Murillo, J. G. Gillmore, *J. Org. Chem.* **2012**, *77*, 6423.

[36] J. M. Guevara-Vela, E. Francisco, T. Rocha-Rinza, A. M. Pendas, *Molecules* **2020**, *25*, 4028.

[37] R. Bader, *Atoms in molecules: a quantum theory*, Oxford University Press, Oxford **1990**.

[38] C. Morris, N. M. Kriege, F. Bause, K. Kersting, P. Mutzel, M. Neumann, *Tudataset: A collection of benchmark datasets for learning with graphs*. arXiv preprint arXiv:2007.08663 **2020**.

[39] B. Ramsundar, P. Eastman, P. Walters, V. Pande, K. Leswing, Z. Wu, *Deep Learning for the Life Sciences*, O'Reilly Media, Sebastopol **2019** https://www.amazon.com/Deep-Learning-Life-Sciences-Microscopy/ dp/1492039837

[40] S. Kearnes, K. McCloskey, M. Berndl, V. Pande, P. Riley, *J. Comput.-Aided Mol. Des.* **2016**, *30*, 595.

[41] Rdkit, Open-source cheminformatics. https://www.rdkit.org

[42] D. Weininger, *J. Chem. Inf. Comput. Sci.* **1988**, *28*, 31.

[43] M. Rupp, A. Tkatchenko, K.-R. Müller, O. A. von Lilienfeld, *Phys. Rev. Lett.* **2012**, *108*, 058301.

[44] M. J. Vainio, M. S. Johnson, *J. Chem. Inf. Model.* **2007**, *47*, 2462.

[45] E. J. Baerends, O. V. Gritsenko, R. van Meer, *Phys. Chem. Chem. Phys.* **2013**, *15*, 16408.

[46] K. P. Murphy, *Machine Learning: A Probabilistic Perspective*, MIT Press, Cambridge, MA **2012**, p. 492.

[47] G. Hoffmann, M. Balcilar, V. Tognetti, P. Héroux, B. Gaüzère, S. Adam, L. Joubert, *J. Comput. Chem.* **2020**, *41*, 2124.

[48] D. Haussler, Convolution kernels on discrete structures. Tech. Rep. (Technical report, Department of Computer Science, University of California at Santa Cruz **1999**.

[49] K. M. Borgwardt, H.-P. Kriegel, *Data Mining, Fifth IEEE International Conference*, IEEE, Washington, DC **2005**, p. 8.

[50] F. Suard, A. Rakotomamonijy, *Esann* **2007**, *355*.

[51] B. Gaüzère, P.-A. Grenier, L. Brun, D. Villemin, *Pattern Recognit.* **2015**, *48*, 356.

[52] N. Shervashidze, P. Schweitzer, E. J. v. Leeuwen, K. Mehlhorn, K. M. Borgwardt, *J Mach Learn Res* **2011**, *12*, 2539.

[53] N. M. Kriege, F. D. Johansson, C. Morris, *Appl Netw Sci* **2020**, *5*, 1.

[54] K. Borgwardt, E. Ghisu, F. Linares-López, L. O'Bray, B. Rieck, Graph kernels: State-of-the-art and future challenges. arXiv preprint arXiv: 2011.03854 **2020**.

[55] S. Ghosh, N. Das, T. Gonçalves, P. Quaresma, M. Kundu, *Comput Sci Rev* **2018**, *27*, 88.

[56] L. Jia, B. Gaüzère, P. Honeine, *Expert Syst Appl* **2022**, *189*, 116095.

[57] B. Schölkopf, A. J. Smola, *Learning with kernels: support vector machines, regularization, optimization, and beyond*, MIT press, Cambridge, MA **2002**.

[58] J. Hartmanis, *Siam Rev* **1982**, *24*, 90.

[59] K. Riesen, H. Bunke, *Image Vision Comput* **2009**, *27*, 950.

[60] L. Jia, B. Gaüzère, F. Yger, P. Honeine, *Structural, Syntactic, and Statistical Pattern Recognition: Joint IAPR International Workshops*, S+SSPR 2020, Padua, Italy, January 21-22, 2021, *Proceedings*, Springer, Berlin **2021**, p. 238.

[61] B. Gaüzère, S. Bougleux, L. Brun, *Joint IAPR International Workshops on Statistical Techniques in Pattern Recognition (SPR) and Structural and Syntactic Pattern Recognition (SSPR)*, Springer, Berlin **2016**, p. 496.

[62] S. Bougleux, B. Gaüzère, L. Brun, 2016 **23rd International Conference on Pattern Recognition (ICPR)**, IEEE Computer Society, Washington, DC **2016**, p. 1701.

[63] D. B. Blumenthal, N. Boria, J. Gamper, S. Bougleux, L. Brun, VLDB **J** **2020**, *29*, 419.

[64] H. Bunke, G. Allermann, *Pattern Recognit Lett* **1983**, *1*, 245.

[65] N. S. Altman, *Am Stat* **1992**, *46*, 175.

[66] M. Fey and J. E. Lenssen, *ICLR Workshop on Representation Learning on Graphs and Manifolds openreview.net* (**2019**).
[67] L. Wu, P. Cui, J. Pei, L. Zhao, X. Guo, Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, Association for Computing Machinery, New York 2022, p. 4840.

[68] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, S. Y. Philip, *IEEE Trans Neural Netw Learn Syst* **2020**, *32*, 4.

[69] J. Zhou, G. Cui, S. Hu, Z. Zhang, C. Yang, Z. Liu, L. Wang, C. Li, M. Sun, *AI Open* **2020**, *1*, 57.

[70] W. L. Hamilton, *Learning* **2020**, *14*, 1.

[71] J. Gilmer, S. S. Schoenholz, P. F. Riley, O. Vinyals, G. E. Dahl, *International conference on machine learning*, PMLR, London **2017**, p. 1263.

[72] T. N. Kipf, M. Welling, Semi-supervised classification with graph convolutional networks. arXiv preprint arXiv:1609.02907 **2016**.

[73] M. Zhang, Z. Cui, M. Neumann, Y. Chen, *Proceedings of the AAAI conference on artificial intelligence*, AAAI Press, Washington, DC Vol. 32 **2018**.

[74] K. Xu, W. Hu, J. Leskovec, S. Jegelka, How powerful are graph neural networks? *arXiv preprint arXiv:1810.00826* **2018**.

[75] P. Velickovic, G. Cucurull, A. Casanova, A. Romero, P. Lio, Y. Bengio, *Stat* **2017**, *1050*, 10.

[76] Y. Shi, Z. Huang, S. Feng, H. Zhong, W. Wang, Y. Sun, Masked label prediction: Unified message passing model for semi-supervised classification. *arXiv preprint arXiv:2009.03509* **2020**.

[77] GCC, Online documentation-GNU Project. https://gcc.gnu.org/onlinedocs/

[78] S. Behnel, R. Bradshaw, C. Citro, L. Dalcin, D. Seljebotn, K. Smith, *Comput Sci Eng* **2011**, *13*, 31.

[79] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, *Adv Neural Inf Process Syst* **2019**, *32*.

[80] Github, The code for these ml experiments is available at the following repository. https://github.com/jajupmochi/RedoxPrediction/

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** L. Jia, É. Brémond, L. Zaida, B. Gaüzère, V. Tognetti, L. Joubert, *J. Comput. Chem.* **2024**, *45*(28), 2383. https://doi.org/10.1002/jcc.27380