GIMADIEV Timur

Modèles prédictifs pour les paramètres cinétiques et thermodynamiques des réactions chimiques

Résumé

Ce travail est consacré à la modélisation QSPR des propriétés cinétiques et thermodynamiques des réactions chimiques à l'aide de l'approche Graphe Condensé de Réaction (CGR). Le CGR permet de coder des structures de réactifs et de produits en un seul graphe moléculaire pour lequel des descripteurs moléculaires peuvent être générés. Une base de données contenant plus de 11000 réactions collectées manuellement a été développée puis utilisée dans la modélisation. Les modèles prédictifs ont été construits pour les constantes de vitesse de réactions Diels-Alder, SN2 et E2 ainsi que pour les constantes d'équilibre des transformations tautomères. Ils sont rendus publics via un portail WEB. Une partie de la thèse concerne une étude de mécanique quantique des réactions entre des sydnones et des alcynes contraints pour lesquels la taille du jeux de données n'était pas suffisante pour produire des modèles statistiquement significatifs.

Résumé en anglais

This work is devoted to QSPR modeling of kinetic and thermodynamic properties of chemical reactions using the Condensed Graph of Reaction (CGR) approach. CGR allows encoding structures of reactants and products into one sole molecular graph for which molecular descriptors can be generated. A comprehensive database containing some 11000 manually collected reactions has been developed then used in the modeling. Predictive models were built for rate constants of Diels-Alder, SN2 and E2 reaction as well as for equilibrium constants of tautomeric transformations. They are available for the users via WEB portal. A part of the thesis concerned quantum mechanics studies of reactions between sydnones and strained alkynes for which the size of the dataset was not sufficient to produce statistically meaningful models. TAU -Tautomeric equilibria RBF -Radial based function BCN -bicyclo-[6.1.0]-nonyne TMTH -3,3,6,6-tetram-ethylthiaheptyne BARACbiarylazacyclooctynone DIBAC -dibenzoazacyclooctyne TCO -trans-cyclooctene DIFO -difluorocyclooctyne Cp(1,3) -1,3-dimethylcyclopropene Cp(3,3) -3,3-dimethylcyclopropene

Table of contents

List of Abbreviations

Résumé en français Introduction

La plupart des approches chémoinformatiques pour l'analyse de données, la visualisation et la modélisation sont développées pour des molécules individuelles codées par des vecteurs de descripteurs. Dans ce contexte, une réaction chimique est un objet complexe, car il implique plusieurs structures moléculaires de deux types -les réactif(s) et les produit(s). Le rendement et les paramètres cinétiques et thermodynamiques de la réaction dépendent des conditions expérimentales (solvant, température, catalyseur, etc…) qui doivent aussi être pris en compte dans la modélisation. Dans ce projet, deux méthodologies différentes sont appliquées pour réduire ces complexités. Afin de réduire la complexité structurelle, on a utilisé l'approche des Graphes Condensés de Réaction (Condensed Graph of Reaction (CGR)), qui combine les structures des réactifs et des produits dans un graphe moléculaire unique, une sorte de pseudo-molécule (Figure 1), pour laquelle des descripteurs moléculaires peuvent être calculés. De plus, quelques descripteurs pertinents pour le solvant et la température ont été proposés afin de réduire la complexité des conditions réactionnelles. Ces méthodologies ont été appliquées ici pour construire des modèles prédictifs des constantes de vitesse pour trois classes de réactions ainsi que pour les constantes d'équilibre de 11 classes de réactions tautomériques.

Malgré tous nos efforts, aucun modèle statistique robuste n'a été obtenu pour les constantes de vitesse des réactions dites "bioorthogonales" entre les sydnones et les iminosydnones avec des cycloalcynes. Par conséquent, la méthode quantique DFT a été appliquée pour accéder aux états de transition des réactions et aux énergies d'activation associées. Toutes les réactions ont été standardisées suivant le processus basé sur la représentation linéaire des CGRs, appelé « signatures CGR ». Au total, 10998 entrées de données brutes de réactions de substitution nucléophilique bimoléculaire (SN2), d'élimination bimoléculaire (E2), de Diels-Alder (DA) et d'équilibre tautomérique (TAU) ont été collectées (Table 1). La colonne « données collectées » dans la Table 1 représente le nombre total de réactions collectées. La colonne « Données standardisées » montre le nombre de réactions après standardisation structurale, correction AAM et suppression des duplicats complets, c'est-à-dire des entrées où tous les champs coincident. Cette dernière raison, étant majeure, est causée par le fait que ce jeu de donées a été collecté simultanément par plusieurs personnes, et parfois, par erreur, une même réaction a pu être extraite deux fois.

Les erreurs dans la structure ou l'abscence de certains champs obligatoires ont été les secondes raisons de la suppression de points de données. La colonne « Jeu modèle » contient des réactions sélectionnées pour la modélisation. La procédure de standardisation de données implique dans ce cas de standardiser, d'identifier les duplicats et de moyenner les propriétés des duplicats. Le nombre total de types de transformation (combinaison de réactifs et de produits) est donné dans le champ « transformation dans le jeu modèle ». La différence entre la taille du jeu modèle et le nombre de transformations reflète le fait que les propriétés ont été mesurées dans plusieurs conditions pour certaines transformations. Modélisation structure-activité des constantes de vitesse des réactions SN2, E2 et Diels-Alder.

Pour chaque réaction, le vecteur de descripteurs résulte de la concaténation des descripteurs caractérisant les structures chimiques, le solvant, et la température.

Les structures chimiques ont été encodées par des descripteurs ISIDA, correspondant à des sous-graphes CGR de différentes tailles et topologies. Chaque solvant a été caractérisé par 15 descripteurs physico-chimiques. La température inverse (1/T) a été aussi utilisée comme descripteur. Pour chaque jeu de données, le flux de travail a impliqué les étapes suivantes : [START_REF] Palm | Tables of Rate and Equilibrium Constants of Heterolytic Organic Reactions[END_REF] préparation des CGRs, [START_REF] Varnek | Substructural fragments: an universal language to encode reactions, molecular and supramolecular structures[END_REF] génération des descripteurs ISIDA [START_REF] Varnek | Substructural fragments: an universal language to encode reactions, molecular and supramolecular structures[END_REF], (3) préparation des vecteurs de descripteurs de réactions comme une combinaison des descripteurs ISIDA, des descripteurs de solvants et de l'inverse de la température, ( 4) construction et validation des modèles individuels en utilisant la méthode de la Régression à Vecteurs de Support [START_REF] Drucker | Support vector regression machines[END_REF], [START_REF] Marenich | Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions[END_REF] sélection des 10 meilleurs modèles individuels et (6) estimation du consensus des 10 modèles au sujet de la propriété modélisée. Dans l'ensemble, 616 types de fragmentation ISIDA différents ont été utilisés; chacun d'eux a conduit à un modèle SVR individuel dont la performance a été évaluée par une validation croisée en 5 paquets répétée 10 fois. Les hyper-paramètres de la SVR et les meilleurs types de fragmentation ont été sélectionnés par un algorithme génétique optimisant le coefficient de détermination estimé en validation croisée. Les performances des modèles consensus de l'estimation des constantes de vitesse (logk) sont données dans la Table 1. On constate que la précision des prédictions de logk est proche de l'erreur expérimentale estimée (0.7-1.0 unités de log).

Une attention particulière a été donnée au jeu de données SN2 qui assemble les réactions effectuées dans 44 solvants différents et impliquant à la fois des nucléophiles anioniques et neutres (voir Table 2). En essayant d'améliorer la performance des prédictions de logk, certains modèles "locaux" correspondant à un solvant particulier ou impliquant un type de nucléophile particulier ont été préparés. Leurs performances étaient néanmoins similaires à celles du modèle "global" (voir Table 2). La visualisation et l'analyse de l'espace de réactions avec l'approche de Cartographie Topographique Générative (« Generative Topographic Mapping », GTM) a été réalisée avec les descripteurs fragmentaux. La puissance des cartes dans l'analyse de données est due à la possibilité de colorer les objets en fonction de différents critères. Ainsi, nous avons coloré les réactions par rapport à leur coeur, aux substrats, et à la nature nucléophile (Figure 2). On voit bien que les cartes séparent bien différentes classes de réactions.

QSPR et modélisation par chimie quantique des constantes d'équilibre tautomérique

Par convention, le logarithme de la constante d'équilibre tautomérique (logKT) est évalué comme la différence entre les valeurs de pKa d'une paire structures tautomères. En prenant en compte l'effet cumulatif des erreurs, ceci pourrait significativement affecter la précision des prédictions de logKT estimés à partir de modèles de pKa. Ici, pour obtenir des modèles pour le logKT, un équilibre tautomérique a été encodé par un Graphe Condensé de Réaction, ce qui nous a permis d'appliquer la procédure de modélisation décrit dans la section 2.2.

Le jeu de données de modélisation contenait 11 types d'équilibres (kéto/énol, amino/imino, azo/hydrazone, pyridol/pyridone, pyridinoid/pyridonoid, phenolimine/keto-amine, thione-enol/keto-thiol, amine-thione/imine-thiol, nitro/acide, forme neutre/zwitterion et cycle/chaîne) mesurés à différentes températures et dans différents solvants. La précision de la prédiction obtenue en validation croisée (RMSE = 0.94 unités de log, Table 2) est similaire aux erreurs expérimentales.

Dans un but de comparaison, des calculs de chimie quantique ont été effectués sur deux jeux de données de test (TEST1 et TEST2) contenant respectivement 23 et 24 équilibres tautomériques. Les valeurs de logKT ont été calculées comme la différence des énergies libres de formation des tautomères estimées par des calculs DFT dans B3LYP/6-311++G(d,p) implémentés dans le logiciel Gaussian09. Les effets de solvant ont été évalués en utilisant le formalisme IEF-PCM du modèle de solvant continu polarisable de Tomasi [START_REF] Tomasi | The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level[END_REF]. Les paramètres SMD [START_REF] Marenich | Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions[END_REF] pour la part non-électrostatique de l'énergie de solvatation ont été utilisés. Nos modèles SVR surpassent significativement ces calculs DFT (Table 3). Il a été aussi montré que nos modèles sont plus performants que l'outil commercialisé par la société ChemAxon, qui fait référence actuellement. Nombre de points de données (Neq), coefficients de détermination (R 2 ), erreurs quadratique moyenne (RMSE en unités de log), et taux de succès de prédiction du tautomère dominant (MT, %).

Cinétique des réactions "bioorthogonales" impliquant des sydnones.

Dans le cadre du projet ANR ClickReal, nous devions construire des modèles pour les constantes de vitesse de réaction (k) des sydnones (SYD) avec des cycloalcynes. Ces réactifs pouvaient potentiellement être utilisés pour des réactions bioorthogonales. Le jeu de données contenait 18 valeurs de k mesurées par nos partenaires. Nous ne sommes pas parvenus à construire un modèle statistique satisfaisant pour logk, quels que soient les descripteurs ou la méthode d'apprentissage machine utilisés. Par conséquent, une série de calculs de DFT en phase gazeuse a été effectuée afin d'identifier les états de transition des réactions et les énergies d'activation associées. L'ensemble de fonctionnelles de Perdew-Burke-Ernzerhof (PBE) et l'ensemble de base 3z implémentés dans le programme Priroda14 ont été utilisés. En accord avec les données expérimentales, le chemin réactionnel contenait 2 états de transition : l'un d'entre eux (TS1, Figure 3) définissait la vitesse de réaction. Les énergies libres d'activation calculées (∆∆G) corrélaient bien avec les observations expérimentales de logk, en utilisant l'équation d'Arrhenius (Figure 4). La précision des prédictions de ∆∆G is 1.97 kcal/mol (RMSE) pour SYD.

Table 4. Reactions of sydnones (A,B,C) with bicyclo-[6.1.0]-nonyne (BCN) and [START_REF] Drucker | Support vector regression machines[END_REF][START_REF] Drucker | Support vector regression machines[END_REF]6, No. R X Cycloalkyne Rate constant (tolerance), M - 1 *sec -1

Reference

A p-Me C6H4 H BCN 0.032 (0.001) [6] B p-Me C6H4 F BCN 42 [6] C p-Me C6H4 F TMTH 1500 [6] Figure 3. Voie réactionnelle pour la réaction A (ligne en pointillés rouges), B (ligne noire pleine), et C (ligne en pointillés bleue), selon la Table 4. Les énergies libres relatives à 298K des molécules par rapport aux réactifs sont montrées. La structure des réactifs, les états de transition et intermédiaires pour la réaction B sont montrés.

Le substituant R du sydnone ainsi que presque tous les atomes du BCN ont été omis pour des soucis de clarté.

la Figure 3 montre que l'étape limitante des réactions est le premier état de transition (TS1), l'intermédiaire étant plus bas en énergie libre que les réactifs de 21 kcal/mol. Selon notre calcul, l'intermédiaire est très instable étant donné qu'il est séparé des produits par un état de transition (TS2) avec une petite barrière (environ 1 kcal/mol). La décomposition intermédiaire est exergonique d'environ 85 kcal/mol. La même chose est vraie pour la réaction beaucoup plus lente A. L'effet du solvant calculé par le modèle IEF-PCM [START_REF] Tomasi | The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level[END_REF] n'a presque aucune influence sur la barrière de réaction, à cause de la compensation. Dans la réaction C entre le fluorosydnone et le TMTH, l'état de transition n'a pas été localisé du tout.

Cette petite barrière d'énergie est pathologique d'une décomposition totalement dominée par les effets entropiques très difficiles à estimer (illustrée Figure 3).

Implémentation des modèles.

Les modèles développés pour les réactions SN2, E2 Chapter 1.

Introduction

Reactions are the main tool for the chemist. New substances cannot be produced without being involved into this process of matter transformation. There are a lot of empirical rules that should be memorized by a chemist in many various examples of transformations. Consequently, this process makes every chemist to a very focused specialist in a narrow area of chemistry. The creation of a chemoinformatic tool helping chemist to extract and store these rules from known experimental data will improve synthesis flexibility of research labs. Chemoinformatics deals with information management and prediction of different properties of various chemical systems [START_REF] Varnek | Chemoinformatics as a Theoretical Chemistry Discipline[END_REF]. A lot of articles are devoted to modeling of different features of molecules, but much less has been written about predicting of properties of reactions like rate constant, conditions, selectivity, yield, etc.

Recently, a gain of interest of chemoinformatics specialists for chemical reaction modeling is observed. It could be related to the fact that, first, a wealth of data in chemical reaction has been accumulated by now and, second, chemoinformatics technologies were developed enough to address new challenges. Several reviews were published [START_REF] Engkvist | Computational prediction of chemical reactions: current status and outlook[END_REF][START_REF] Baskin | Artificial intelligence in synthetic chemistry: achievements and prospects[END_REF]. It was stated that organic synthesis is a rate-limiting factor in drug discovery and the usage of artificial intelligence tools could revolutionize medicinal chemistry [START_REF] Blakemore | Organic synthesis provides opportunities to transform drug discovery[END_REF].

The development of chemoinformatics approaches for reactions is limited by the fact that it is a complex process: it involves reagent(s) and product(s), depends on their concentrations, reaction conditions (like temperature, pressure, radiation), This makes representation of reactions in electronic databases much more problematic than the storage of simple molecules, requiring specific procedures for standardization. Another rising problem is incompleteness of experimental data. There is a lot of information about reactions, but popular databases like Reaxys and SciFinder provide very heterogeneous information. Often, only one major component of reaction is reported. The information about kinetics of reactions is rare and usually comes from only one source, so there is no cross-validation of values between various experimental methods.

Currently, the main tool of reactions studies is quantum chemistry. Chapter 2.

Review on chemical reaction modeling studies

Thermodynamics and kinetics of reaction

A chemical equation is a special symbolic representation of a chemical reaction.

The reactant entities are given on the left-hand side and the product entities on the righthand side.

𝑥𝐴 + 𝑦𝐵 + ⋯ = 𝑛𝐶 + 𝑚𝐷 + ⋯

The numbers next to the formulae of entities represent number of molecules participating reaction act are called stoichiometric coefficients [START_REF] Cohen | Quantities, Units and Symbols in Physical Chemistry[END_REF]. The first-ever chemical equation was diagrammed by Jean Beguin in 1610 in first edition of "Tyrocinium Chymicum" book [START_REF] Beguin | Tyrocinium Chymicum[END_REF].

The reaction is a complex instance that requires a description of physical conditions and chemical environment of compounds along with description of reagents and products. Chemical environment includes compounds such as solvents, catalysts, catalytic poisons, additives, etc, which do not participate in reaction according to its chemical equation, but can influence it. Physical conditions include physical factors that are significant for reaction: temperature, pressure, irradiation, etc. All of them influence the macroscopic parameters characterizing the reaction process : yield, rate, selectivity, equilibrium constant etc. From chemists' point of view, the most important endpoint of reactions is the yield of desired compound. The latter depends on the speed of product and byproducts accumulation (reaction rate) and the equilibrium constant, the key thermodynamic parameter of a chemical process.

Properties of reactions can be divided into thermodynamic and kinetic ones.

Thermodynamic properties reflect the position of equilibrium, heat of reaction and work against external forces. Kinetic properties reflect speed of reagents conversion or product accumulation.

Among the most important thermodynamic properties of chemical reaction are enthalpy, entropy and free energy of reaction. The enthalpy, H, comprises a system's internal energy, which is the energy required to create the system, plus the amount of work required to make room for it by displacing its environment and establishing its volume and pressure. The difference between enthalpy of products and reagents is equal to the heat released or absorbed in reaction, provided that system has constant pressure:

Q = -ΔH.
Entropy, S, is a thermodynamic function reflecting the number of microscopic configurations that a thermodynamic system can have in a state with defined macroscopic variables. In other words, it measures the degree of disorder in the thermodynamic system. The difference between enthalpy (total energy) and entropy (amount of useless energy accounting for disordered movement) defined Gibb's free energy G = H -TS, reflecting the amount of work that a thermodynamic system can perform (where T is absolute temperature of the system, in Kelvin). The standard change of Gibb's free energy in a reaction (Eq. 1), under isothermal and isobaric conditions is related with the equilibrium constant, equation (Eq. 2), K, as shown in equation (Eq. 3).

Eq. 1 ∆𝐺 = ∆𝐻 -𝑇∆𝑆

Eq. 2 𝐾 = [456789:;][456789:=]… [?@AB@C:;][?@AB@C:=]… Eq. 3 ∆𝐺 = -𝑅𝑇 • log 𝐾
where R is the gas constant.

The equilibrium constant defines the relative concentrations of products and reagents at chemical equilibrium. Thus, equation ( 2) and ( 3) used for calculation of reagent and product concentration. While direct experimental measurement of ΔG is rather challenging, equilibrium constant could be readily calculated from (2) if composition of mixture at equilibrium is known.

Chemical kinetics deals with the study of reaction evolution in time. Eq. 4

𝑘 = 𝐴𝑒 K L M NO
where A is pre-exponential factor, and E a is activation energy of reaction.

Thermodynamics and kinetics are closely connected with each other, e.g. rate constants of forward, kf, and backward reaction, kb, can be used to calculate equilibrium constant for this process, K= kf / kb.

Competitive reversible reactions could proceed under thermodynamic or kinetic control. A reaction selectivity will highly depend on the dominating type of control that itself depends on reaction condition. Here (Figure 6), product P1 has lower activation energy for reaction of its formation (greater reaction rate), but P2 is more stable (dominate at equilibrium). In this case in the beginning of reaction the product P1 is formed more rapidly but over time the amount of P2 rises and finally it becomes dominating.

Figure 6. Thermodynamic and kinetic reaction control.

For example, high temperature and catalysts shorten the time required for reaching equilibrium and favor the thermodynamically controlled product. With properly selected reaction time, the kinetically controlled product can dominate in mild conditions.

Knowing reaction rate constants and equilibrium constant or ΔG of reaction, one could shift reaction in desired direction and maximize the yield of desired product.

Unfortunately, the measure of kinetic and thermodynamic data for reactions is quite expensive, and prediction tools are either too computationally demanding (like quantum chemical prediction), or insufficiently developed (like chemoinformatics approaches), or have limited applicability (like LFER approaches).

Prediction of reaction characteristics for synthesis optimization is an unsolved problem. Quantum chemical methods could hardly be used as a priori calculation tools mainly due to their time-and resource-consummation. However, in 2013 scientists proposed an approach for selection of the best solvent in bimolecular substitution reaction, combining quantum mechanical computations of the reaction rate constant in a few solvents with linear regression based on solvent descriptors [START_REF] Struebing | Computer-aided molecular design of solvents for accelerated reaction kinetics[END_REF]. According to authors, it took several days to make one prediction.

Although properties of pure compounds are widely predicted by means of Quantitative Structure-Property (Activity) Relationship (QSPR, QSAR) approach [START_REF] Palm | Tables of Rate and Equilibrium Constants of Heterolytic Organic Reactions[END_REF] that favorably differs from quantum chemical calculations in terms of required effort, only few works were devoted to QSPR modeling of reactions. The main reason is that reaction is a more complex ITEM for modeling than a molecule [START_REF] Baskin | Artificial intelligence in synthetic chemistry: achievements and prospects[END_REF]. First, a reaction involves several molecules. Next, a reaction has a sense, which must be accounted for in predictions: the backward process from products to reagents is characterized by DGback=-DG and Kback=1/K. Moreover, some chemical compounds such as solvents, catalysts and additives influence the reaction but they are usually omitted from the reaction equation.

And finally, physical conditions (temperature, pressure) are very important for some properties. Hence, account for all this effect is required in QSAR/QSPR study of reactions and special approaches for handling reaction complexity need to be introduced.

In spite of the huge amount of known reactions, only a tiny portion has any kinetic or thermodynamic data and this information has not been transferred from primary sources to electronic databases. There is no publicly available database of kinetic characteristics of diverse reactions in common conditions (there are some small databases of gas phase reactions [START_REF] Manion | NIST Chemical Kinetics Database[END_REF][START_REF]Mellouki W Chemical Kinetics Database on oxygenated VOCs gas-phase reactions[END_REF]). The largest databases like Reaxys and CAS REACT annotate only the fact that kinetics of thermodynamic data were measured but not the values.

Reaction representation

As already mentioned, reaction description has to include information about reagents (structural information), and products, temperature, solvent, catalyst,

respectively (condition descriptors).

There are several ways to represent structural information of reactions, based on graph representation (usually represented as connection tables), linear notations systems (reaction is represented by an alphanumerical array) or descriptor representation calculated usually on the basis of previous two. The numeric encoding of its characteristics under the form of a vector of descriptors is key feature of chemoinformatics.

Generally, one can define 3 main ways to represent a reaction in chemoinformatics:

(i) as ordered set of reagent and product molecules, (ii) reaction center based representation and (iii) product-reagent difference. In our work we proposed the special type of reaction representation [START_REF] Polishchuk | Structure-reactivity modeling using mixture-based representation of chemical reactions[END_REF] in which reagent and product part of reaction is treated as mixture with given composition [START_REF] Oprisiu | QSPR Approach to Predict Nonadditive Properties of Mixtures. Application to Bubble Point Temperatures of Binary Mixtures of Liquids[END_REF].

Reagent-product representation

Reagent-product representation of chemical reaction is based on consideration of common reaction equation: reagents in reaction are enumerated and followed by products. Using special approaches reagents are explicitly separated from products to avoid confusion.

One of the most popular system for molecule representation is the simplified molecular-input line-entry system (SMILES), developed by D. Weininger [START_REF] Weininger | SMILES, a chemical language and information system. 1. Introduction to methodology and encoding rules[END_REF]. It is an ASCII string that is compiled according to specification that encodes the structure of chemical species. The main advantage of this type of representation is saving of space for storage, machine and human readability. An algorithm of canonicalization [START_REF] Weininger | SMILES. 2. Algorithm for generation of unique SMILES notation[END_REF] is used for generation of a unique SMILES representation of molecule. As a development of SMILES, the SMIRKS [START_REF]SMIRKS -A Reaction Transform Language[END_REF] represents reaction as transformation of reagents SMILES into products SMILES. Reagents and product molecules are separated by dots as usually for disconnected components in SMILES. After list of reagents ">" symbol opens list of additives which follows another ">" symbol, after which products are specified. In the absence of additives, the ">>" symbol separates the lists of reagents and products. Atoms and their correspondence (atom-to-atom mapping, see below) could be specified using numbers (see Figure 7). Chemical table file (CT File) formats developed by Molecular Design Limited [START_REF] Dalby | Description of several chemical structure file formats used by computer programs developed at Molecular Design Limited[END_REF] represent the most popular file formats for storage and exchange by structural information. For reaction representation, the RXN file format was developed. It is based on enumeration of specified number of reagents followed by enumeration of products.

Molecules are represented in MOL format. RDF file format is an extension of RXN providing functionality to store factual data associated with a given reaction in text, numeric and structure (for reagents, catalysts, etc) fields. Other popular file format based on Chemical Markup Language has an extension for reaction representation [START_REF] Holliday | Chemical Markup, XML, and the World Wide Web. 6. CMLReact, an XML Vocabulary for Chemical Reactions[END_REF].

Reaction center representation

Reagent-product representation is not an effective way for reaction storage, since every atom is indeed represented twice -in reagents and in products. Reaction center type representations avoid this drawback, encoding the changes occurring in reaction.

This type of representation is tightly connected with the notion of reaction center, i.e. the atoms incident to edges associated with bond order changes). With a few exceptions, changes in reaction are mostly associated with bond formation, cleavage and bond order changes.

Detection of reaction center is a crucial step required for this type of representations. Atom-to-atom mapping (AAM) or bond-to-bond mapping (BBM) is the essential procedure used for identification of reaction center [START_REF] Chen | Automatic reaction mapping and reaction center detection[END_REF]. AAM/BBM establishes one-to-one correspondence between atoms/bonds of reagents and products assigning each atom a unique label (Figure 9), AAM/BBM allows identification of atoms with altered environment and thus enables automated reaction center detection. For reaction on Figure 9, three atoms comprise reaction center Br 1 , C 2 , N 9 . Reaction center can be annotated as the simplest possible reaction of a given type: Br-C + N = CH-N + Br - (implicit hydrogens are omitted). Automated, computer-based AAM is one of the key problems in reactions processing. However the procedure itself is error prone and rather computationally intensive [START_REF] Chen | Automatic reaction mapping and reaction center detection[END_REF]. Numbers represent atom-to-atom mapping.

Several reaction center-based representations were proposed. Initially they were proposed for reaction classification rather than for information storage and retrieval.

Reaction centers are indeed the "signatures" that differentiate reaction [START_REF] Chen | Reaction Classification and Knowledge Acquisition[END_REF]. They were used for creation of an ontology of chemical reactions and creation of unique reaction identifiers, like the IUPAC nomenclature for chemical compounds. The first reaction classification scheme based on electron redistribution in pericyclic reactions including 6 atoms in reaction center was proposed by Balaban [START_REF] Balaban | Chemical graphs. 3. Reactions with cyclic 6-membered transition states[END_REF]. Arens proposed an approach for textual representation of bond types valid both for cyclic (when reaction center forms cyclic graph) and linear reaction centers (without cycles) [START_REF] Arens | A formalism for the classification and design of organic reactions. I. The class of (-+)n reactions[END_REF][START_REF] Arens | A formalism for the classification and design of organic reactions. II. The classes of (+ -)n + and (-+)n -reactions[END_REF][START_REF] Arens | A formalism for the classification and design of organic reactions III. The class of (+ -)nC reactions[END_REF], Figure 10. Hendrickson developed Balaban's approach and created coherent classification system of pericyclic reactions based on bonds redistribution within 4-, 5-and 6-membered rings [START_REF] Hendrickson | The Variety of Thermal Pericyclic Reactions[END_REF]. The latter was further extended to other reaction types and called "comprehensive system for classification and nomenclature of organic reactions" [START_REF] Hendrickson | Comprehensive System for Classification and Nomenclature of Organic Reactions[END_REF], Figure 10. Zefirov and Tratch developed a hierarchical system for reaction representation and a classification that considers not only the reaction center, but its higher level of generalization [START_REF] Zefirov | An Approach to Systematization and Design of Organic Reactions[END_REF][START_REF] Tratch | A Hierarchical Classification Scheme for Chemical Reactions[END_REF].

Albeit universal and widely accepted methodologies of reaction classification haven't been developed, the efforts were not made in vain and finally led to progress in reaction data mining. One of the most promising approaches in this emerging field of application is the Condensed Graph of Reaction methodology. The Condensed Graph of Reaction approach is based on the early work of Vladutz [START_REF] Vléduts | Concerning one system of classification and codification of organic reactions[END_REF], who used it for reaction information storage and classification. He was the first to give the definition of reaction center as group of atoms in which changes of bonds take place during the reaction and proposed a method for representing the reaction center in terms of bond changes [START_REF] Vladutz | Do we still need a classification of reactions?[END_REF]. He proposed to merge reaction center atoms from reactant and product graphs, and specially mark changes in the bond orders, (Figure 10). He called the obtained graph the "skeletal scheme of reaction". The latter is common to reactions of the same type and encodes the changes taking place in the reaction center.

Kiho proposed to construct these graph notations not only for the reaction center, but also for the whole reaction [START_REF] Kikho | Formal definition of some notions of quantitative oraganic chemistry[END_REF], however his paper was unnoticed and the same approach was proposed independently by S. Fujita [40]. He called the superimposed graph of reagents and products Imaginary Transition State [40]. The latter had three types of bonds: in-bonds, out-bonds and par-bonds standing for the bonds that are formed, cleaved and unaffected in the reaction, respectively. Thus the whole reaction is represented by a single graph with different labels on the edges that trace bond transformation. Further development of this approach was made in the group of Kauffman, who renamed Imaginary Transition State representation as Condensed Graph of Reaction and distinguished only two types of bonds: ordinary bonds (single, double, etc.) and dynamic bonds (single broken, single created, etc.) [START_REF] Jauffret | Machine learning of generic reactions: 2. toward an advanced computer representation of chemical reactions[END_REF]. To additionally represent changes in stereochemistry and atom formal charges in reaction he added pseudo-atoms standing for atomic charges, stereochemistry, etc. Thus any changes in the reaction were encoded by bond reorganization.

Any atom-to-atom-mapped reagent-product representation of a reaction can be easily converted into CGR, Figure 11. A CGR can be easily represented as any molecular 2D-sketch. The bonds, which are created, broken, or modified during the reaction are called dynamic bonds and represented by specific notations denoting the change, e.g. a single bond broken by the reaction, Figure 11. The bonds that are not changed in reaction (single, double, aromatic, etc.) are represented conventionally. Recently we proposed discontinue support for pseudo-atoms in CGR and instead introduce "dynamic atoms" to represent changes in atomic properties [START_REF] Gimadiev | Assessment of tautomer distribution using the condensed reaction graph approach[END_REF].

Representation based on difference in structure of reagents and products.

While detection of reaction center requires AAM, the changes in reaction could be detected subtracting features or reagents from products. This idea is used in difference reaction representation.

Ugi and Dugundji [START_REF] Dugundji | An algebraic model of constitutional chemistry as a basis for chemical computer programs[END_REF] proposed to represent reaction as difference of bondelectron matrices of products (E matrix) and reagents (B matrix). B and E-matrices of reagents and products, correspondingly, are symmetric matrices where elements of main diagonal represent number of valence electrons on the lone pairs of certain atoms and offdiagonal elements are bond orders (Figure 12). A reaction is represented as a new difference matrix R = E -B. R gives a description of changes taking place in the reaction center and has some other interesting mathematical properties [START_REF] Dugundji | An algebraic model of constitutional chemistry as a basis for chemical computer programs[END_REF]. The Ugi -Dugundji formalism was used to plan organic synthesis and predict reaction pathways [START_REF] Gasteiger | EROS A computer program for generating sequences of reactions[END_REF] and to search for the shortest distance between the reactant and product [START_REF] Jochum | The Principle of Minimum Chemical Distance(PMCD)[END_REF]. AAM is required for R-matrix creation. This type of representation is useful for many applications but has not been used for reaction characteristics modeling.

Reaction descriptors

Graph representations of chemical objects are mainly used for primary structural information storage but they need further preprocessing in order to encode information leading the way to numerical descriptors that are perfect for in silico processing.

Representation of chemical items as vectors of numerical descriptors allows for fast information retrieval and can be used in conjunction with most machine learning methods i.e. QSAR/QSPR modeling. However, descriptor representation is not always invertible, and the structure of the item usually may not be always be restored from descriptor values -which implies that the encoding process may trigger some loss of information.

Descriptors are values that otherwise represent features of a chemical item.

Evolution of QSPR/QSAR modeling resulted in numerous descriptors of different types.

The comprehensive compilation of descriptors were made by Todescini and Consonni [START_REF] Todeschini | Handbook of Molecular Descriptors[END_REF]. They can be categorized into several types. For example, according to origin, descriptors could be classified as:

• Measured physicochemical properties such as "logP(o/w)," the logarithm of the octanol/water partition coefficient, molecular polarizability and refraction.

• Calculated descriptors are mathematical models of some kind, from the simplest counting of carbon atoms and summing up of the molecular mass, to descriptors representing predicted properties according to QSPR equations, themselves relying on simpler descriptors.

According to dimensionality of the molecular representation required for calculation of descriptors one can distinguish:

• Оne-dimensional (1D) descriptors that are calculated from the composition of compounds,

• 2D descriptors that are calculated from a planar graph representation of a molecule,

• 3D descriptors such as "molecular volume" or "surface area" that are derived from molecular conformation, and

• 4D descriptors that require ensemble of conformations or molecular dynamics trajectories.

According to application typically the following descriptor types are distinguished:

• Fragment descriptors that are usually represented by binary and integer vectors monitoring the existence or the frequency of fragment occurrence in a structure [START_REF] Baskin | Chapter 1. Fragment Descriptors in SAR/QSAR/QSPR Studies, Molecular Similarity Analysis and in Virtual Screening[END_REF]. The main advantage of fragment descriptors is their universality [START_REF] Baskin | On the Basis of Invariants of Labeled Molecular Graphs[END_REF].

• Topological indices or connectivity indices are planar molecular graph invariants used for numerical representation of graph topology [START_REF] Gonzalez-Diaz | Medicinal Chemistry and Bioinformatics -Current Trends in Drugs Discovery with Networks Topological Indices[END_REF].

• Physicochemical descriptors can be obtained from experimental measurements of compounds physicochemical properties. The most frequently used descriptor of this type is the logarithm of the octanol/water partition coefficient, a measure of the hydrophobic character of a molecule. Nowadays, these characteristics are usually calculated using QSPR models.

• Quantum-chemical descriptors are characteristics that obtained from approximate solution of Schrodinger's equation for molecules [START_REF] Karelson | Quantum-Chemical Descriptors in QSAR/QSPR Studies[END_REF]. Energetic, molecular orbital and electron distribution descriptors are the most widely used.

• Descriptors of molecular fields are the values that approximate interaction of the molecule with some virtual probe. This descriptors are often used in 3D-QSAR modeling of biological activity, for example, in CoMFA approach [START_REF] Cramer | Comparative molecular field analysis (CoMFA). 1. Effect of shape on binding of steroids to carrier proteins[END_REF].

• Pharmacophore descriptors show occurrence of pharmacophore-labelled fragments (usually pharmacophore pairs or triplets with defined distance between centers) in a molecule [START_REF] Bonachéra | Fuzzy Tricentric Pharmacophore Fingerprints. 1. Topological Fuzzy Pharmacophore Triplets and Adapted Molecular Similarity Scoring Schemes[END_REF]. A pharmacophore is the ensemble of steric, electronic and other physico-chemical properties that are necessary to ensure optimal supramolecular interactions with a specific biological target structure.

• Substituent constants, first introduced by Hammet [START_REF] Hammett | The Effect of Structure upon the Reactions of Organic Compounds. Benzene Derivatives[END_REF], reflect electronic or sterical influence of a given substituent on the core molecule.

• Descriptors of molecular similarity report the molecular similarity with respect to some common set of reference compounds.

The descriptors best suited to encode the structural information of reactions are fragment descriptors of reagents/products/CGR. Since fragment descriptors represent per se a vast chapter of possible fragmentation schemes and strategies to capture specific chemical information, this category was the only one exploited in this work (see below).

Quantum-chemical descriptors and substituent constants are often used in conjunction with linear regression to model some reaction properties (usually kinetics) or to reveal reaction mechanism.

Fragment Descriptors

An important advantage of fragment descriptors is related to the simplicity of their calculation, storage and interpretation [START_REF] Varnek | Substructural fragments: an universal language to encode reactions, molecular and supramolecular structures[END_REF][START_REF] Zefirov | Fragmental Approach in QSPR[END_REF][START_REF] Artemenko | Artificial neural network and fragmental approach in prediction of physicochemical properties of organic compounds[END_REF]). They belong to information-based descriptors [START_REF] Jelfs | Estimation of p K a for Druglike Compounds Using Semiempirical and Information-Based Descriptors[END_REF], which encode the information stored in molecular structures. This contrasts with knowledge-based (or semi-empirical) descriptors derived from consideration of the mechanism of action. Owing to their versatility, fragment descriptors can efficiently be used to build structure-property models, perform similarity searches, virtual screening and in silico design of chemical compounds with desired properties.

The following types of fragment descriptors can be distinguished [START_REF] Varnek | Chemoinformatics Approaches to Virtual Screening[END_REF]: Simple Fixed

Types Fragments [START_REF] Smolenski | [END_REF], WLN and SMILES Fragments [START_REF] Wiswesser | How the WLN began in 1949 and how it might be in 1999[END_REF], Sequences [START_REF] Varnek | Substructural fragments: an universal language to encode reactions, molecular and supramolecular structures[END_REF], Atomcentered Fragments [START_REF] Rosenkranz | Mutagens, carcinogens, and computers[END_REF][START_REF] Benson | Additivity Rules for the Estimation of Molecular Properties. Thermodynamic Properties[END_REF], [START_REF] Zhang | Structure-Based Classification of Chemical Reactions without Assignment of Reaction Centers[END_REF], Atom Pairs and Topological Multiplets [START_REF] Bonachéra | Fuzzy Tricentric Pharmacophore Fingerprints. 1. Topological Fuzzy Pharmacophore Triplets and Adapted Molecular Similarity Scoring Schemes[END_REF][START_REF] Avidon | Structureactivity relationship oriented languages for chemical structure representation[END_REF][START_REF] Schneider | Scaffold-Hopping" by Topological Pharmacophore Search: A Contribution to Virtual Screening[END_REF], Substituents and Molecular Frameworks [START_REF] Free | A Mathematical Contribution to Structure-Activity Studies[END_REF][START_REF] Hansch | p -σ-π Analysis. A Method for the Correlation of Biological Activity and Chemical Structure[END_REF][START_REF] Fujita | A New Substituent Constant, p, Derived from Partition Coefficients[END_REF],

Basic Subgraphs [START_REF] Randic | Representation of molecular graphs by basic graphs[END_REF], Mined Subgraphs [START_REF] Kramer | Molecular feature mining in HIV data[END_REF][START_REF] Asai | Efficient substructure discovery from large semi-structured data[END_REF], Random Subgraphs [START_REF] Graham | Information Content in Organic Molecules: Quantification and Statistical Structure via Brownian Processing[END_REF], Library Subgraphs [START_REF] Sanderson | Computer Prediction of Possible Toxic Action from Chemical Structure; The DEREK System[END_REF]. Despite many different type of fragment descriptors being proposed, only some of them became widely used. Below, two most universal approaches for fragment descriptor generation are discussed: ISIDA and SiRMS fragments.

ISIDA descriptors

ISIDA (shorthand for In SIlico Data Analysis) fragments -simple way for encoding molecules and reactions. For both molecules and CGR, the ISIDA Fragmentor [START_REF] Varnek | ISIDA -Platform for Virtual Screening Based on Fragment and Pharmacophoric Descriptors[END_REF] produces a vector of integers counting the occurrences of molecular fragments of different topology. There are two types of the ISIDA descriptors: sequences of atoms and/or bonds, and augmented atoms (circular fragments) representing a selected atom with its closest environment, Figure 13. There are several main parameters controlling ISIDA fragment descriptor generation:

• Fragmentation scheme. Three fragmentation schemes are supported: sequences, augmented atoms (circular fragments) and topological triplets. User can specify whether to take atom and bond labels in the generated fragments into account.

For example, in Figure 13 one can see values for sequence descriptors.

• Minimal and maximal length of sequence or augmented atom radius. All sequences having specified length or augmented atoms having given topological distance from central atom are generated. This parameter influences the balance between long and short fragments. If short fragments are too general such description will be non-informative. Long fragments could be too selective and numerous and thus could introduce noise in the model.

• Do_all_ways -flag switches off default calculation of shortest path sequences and fragments generated by all possible detours on the graph. It is relevant only for structures with cycles.

• Formal charge option adds the information of the formal charge on an atom into fragment description. This option is useful to differentiate protonation states of molecules. In case of CGR ISIDA fragment descriptors are universal and were used for many different QSAR/QSPR tasks: prediction of H-bond stability [START_REF] Glavatskikh | Predictive Models for the Free Energy of Hydrogen Bonded Complexes with Single and Cooperative Hydrogen Bonds[END_REF][START_REF] Fiorella | Individual Hydrogen-Bond Strength QSPR Modelling with ISIDA Local Descriptors: a Step Towards Polyfunctional Molecules[END_REF], halogen bonding [START_REF] Glavatskikh | Predictive Models for Halogen-bond Basicity of Binding Sites of Polyfunctional Molecules[END_REF], metal complexation [START_REF] Solov'ev | QSPR ensemble modelling of the 1:1 and 1:2 complexation of Co2+, Ni2+, and Cu2+ with organic ligands: relationships between stability constants[END_REF], and many other biological and physico-chemical properties of molecules [START_REF] Varnek | Substructural fragments: an universal language to encode reactions, molecular and supramolecular structures[END_REF][START_REF] Varnek | ISIDA -Platform for Virtual Screening Based on Fragment and Pharmacophoric Descriptors[END_REF][START_REF] Ruggiu | ISIDA Property-labelled fragment descriptors[END_REF], if to name a few.

SiRMS descriptors

SiRMS (shorthand for SImplex Representation of Molecular Structure) are topological multiplets (see 2.3.1). The descriptor vector encodes frequency of presence of all possible atomic multiplets (atom combinations) having given number of vertices.

MultipletSSS could contain atoms that are not bound, but unlike other topological multiplets the topological distance between atoms is not specified in SiRMS, Figure 14.

Initially, SiRMS contained only so called simplexes (tetratomic multiplets) but later fragments with any number of centers (from two to six, as an option) became used [START_REF] Kuz'min | Hierarchical QSAR technology based on the Simplex representation of molecular structure[END_REF].

The latter are still called simplexes due to historical reasons. Description of simplexes could contain bond orders, stereochemical configuration of simplex, atomic property labels (like charge, lipophilicity, polarizability, etc.). Optionally, simplexes containing atoms not bound to other atoms of a given fragment could be discarded. SiRMS descriptors were shown to be generally applicable for prediction of many different properties: biological activity of compounds [START_REF] Kuz'min | Quantitative structure-affinity relationship of 5-HT1A receptor ligands by the classification tree method[END_REF][START_REF] Muratov | Per aspera ad astra: application of Simplex QSAR approach in antiviral research[END_REF], environmental toxicity [START_REF] Polishchuk | Application of Random Forest Approach to QSAR Prediction of Aquatic Toxicity[END_REF], solubility [START_REF] Kovdienko | Application of Random Forest and Multiple Linear Regression Techniques to QSPR Prediction of an Aqueous Solubility for Military Compounds[END_REF] and other physico-chemical properties [START_REF] Mokshyna | Predictive QSPR Modelling for the Second Virial Coefficient of the Pure Organic Compounds[END_REF]. One of their specific features is the possibility to correctly represent stereochemistry-related properties [START_REF] Polishchuk | Structural, Physicochemical and Stereochemical Interpretation of QSAR Models Based on Simplex Representation of Molecular Structure BT -Advances in QSAR Modeling: Applications in Pharmaceutical, Chemical, Food, Agricultural and Environmental Sciences[END_REF].

Quite recently, mixture SiMRS descriptors [START_REF] Oprisiu | QSPR Approach to Predict Nonadditive Properties of Mixtures. Application to Bubble Point Temperatures of Binary Mixtures of Liquids[END_REF] incorporating calculation of simplexes for chemically disconnected molecules were proposed. For a mixture, one fragment could contain parts coming from several molecules. This approach was successfully used for prediction of properties of binary mixtures [START_REF] Mokshyna | QSPR-Modeling for the Second Virial Cross-Coefficients of Binary Organic Mixtures[END_REF]. The "Quasi-mixture" approach, in which every compound is represented as a mixture with itself, was shown to be superior over the single molecule approach in physical properties prediction [START_REF] Mokshyna | Quasi-Mixture' Descriptors for QSPR Analysis of Molecular Macroscopic Properties. The Critical Properties of Organic Compounds[END_REF].

Particular types of reaction descriptors

The reaction descriptor vector should encode two essential features of chemical reaction -the transformation (i.e. structure of reagents and products) and reaction conditions. To simplify the situation, modeling usually addresses either a set of reactions under common conditions, or a single reaction under variable conditions. In this case, the constant part of descriptor vector is neglected.

In this part we will describe some techniques for reaction descriptor calculation.

Reagent-product descriptors

The simplest way to describe reaction is to consider it as a combination of reagents and products. Reagents and products in this approach are considered as molecular graphs. Thus this approach is based on molecular descriptors making it rather general since it does not require special reaction-oriented descriptors.

Reagent descriptors

Let us consider the case of a series of reactions of the same type. In this case, the structure of product is fully determined by reagents and thus only reagents could be considered.

The Linear Free Energy Approach (LFER) is based on this assumption. Moreover, within LFER, reagents are congeneric, i.e. they differ by substituent. Thus changes in predicted value could be explained only by electronic and steric influence of the substituents. The first substituent descriptor was introduced by Hammett [START_REF] Hammett | The Effect of Structure upon the Reactions of Organic Compounds. Benzene Derivatives[END_REF] and was used for prediction of acidity. Later, Taft [START_REF] Taft | Polar and Steric Substituent Constants for Aliphatic and o-Benzoate Groups from Rates of Esterification and Hydrolysis of Esters1[END_REF] introduced the steric factor of substituents.

Later, a lot of different substituent constants were proposed to describe inductive, mesomeric and other electronic and steric effects. Substituent constants [START_REF] Reichardt | Solvents and Solvent Effects in Organic Chemistry[END_REF] descriptors were used for prediction of reaction rate or equilibrium constants k using linear equation of the type

P P Q = σρ
, where k0 is reaction rate with hydrogen as substituent, s is substituent constant and r is reaction and scaffold type dependent constant. Often, multiparametric LFER equations were built to describe dependence of property under study on several substituent constants or solvent parameters. Since these descriptors have clear physical meaning the dependencies are used for interpretation and to study reaction mechanism. Although the usage of LFER approach for predictive purpose has mostly historical importance, it is still used for mechanism elucidation and interpretation [START_REF] Reichardt | Solvents and Solvent Effects in Organic Chemistry[END_REF][START_REF] Salin | Solvent Effect on Kinetics and Mechanism of the Phospha-Michael Reaction of Tertiary Phosphines with Unsaturated Carboxylic Acids[END_REF] Usage of the approach with predictive purpose is limited since it has a very narrow applicability mainly due to two facts: (i) only congeneric reactions could be predicted and

(ii) descriptors used are obtained experimentally. Nevertheless, from time to time predictive models based on LFER approach are proposed. For example, in the work [START_REF] Struebing | Computer-aided molecular design of solvents for accelerated reaction kinetics[END_REF] LFER equation was used to correlate quantum-chemically predicted activation energies for a given reaction in different solvents with descriptors of solvents to predict optimal conditions for the reaction.

The idea that reagent descriptors could be used to predict reaction properties is sometimes exploited in the field of QSRR. Various descriptors [START_REF] Sukhachev | QSPR approach to the calculation of rate constants of homolysis of nitro compounds in different states of aggregation[END_REF] of compounds including topological indices, information indices based on charge distribution in molecules, fragment descriptors and others were used for prediction of rate of homolysis of nitrocompounds. In this case, focusing the description on anything else but the reagents is impossible, since the mechanism of reaction is unknown and reaction results in complex mixture of products.

In the work of Marcou et al different descriptor creation strategies were benchmarked [START_REF] Marcou | Expert System for Predicting Reaction Conditions: The Michael Reaction Case[END_REF] in order to predict Michael reaction feasibility under given conditions (solvent type and catalyst type). They compared reagent based descriptors (420 CDK descriptors, MOLMAP [START_REF] Zhang | Structure-Based Classification of Chemical Reactions without Assignment of Reaction Centers[END_REF][START_REF] Latino | Genome-scale classification of metabolic reactions and assignment of EC numbers with self-organizing maps[END_REF], ISIDA fragments and EED descriptors [START_REF] Elhabiri | Electrochemical Properties of Substituted 2-Methyl-1,4-Naphthoquinones: Redox Behavior Predictions[END_REF]) and came to conclusion that the models employing them have the same predictive performance that the one based on descriptors with explicit reaction center encoding. They came to conclusion that within the considered training set, condition specific structural "patterns"

or "signatures" can be established even in absence of explicit knowledge of the reaction center itself.

Concatenated descriptors of reagents and products

In case when for a given reagent several products can be formed, and modeling reaction property depends on the product formed, the reagent based approach described above can no longer be used. In this case, the reaction descriptor should explicitly encode both reagent and product features. The natural idea to solve this problem is to concatenate reagent and product descriptor vectors. This approach is quite universal and can be used for any chemical reaction, but leads to rather long descriptor vectors. Large descriptor vectors could introduce noise and spurious correlations into a model. If several reagents and products are involved, they should be represented in some specific order (reaction descriptor would be concatenation of descriptors for A, B, C and D molecule).

Thus, careful curation of the dataset is required.

Kravtsov et al [START_REF] Kravtsov | Prediction of Rate Constants of SN2 Reactions by the Multicomponent QSPR Method[END_REF] used concatenated representation of reagents and products based on topological, physicochemical and quantum chemical descriptors for the modeling reaction rate of SN2 reaction. The transformation descriptor vector was concatenated with vector of solvent descriptors and reaction temperature. Using descriptor selection on the basis of Fast Stepwise Multiple Linear Regression and artificial neural networks as machine learning method the first chemoinformatics model able to predict reaction rate in different conditions were built. Similar approach was used in the work of the same authors for classification of preferable mechanism of nucleophilic substitution (SN1 or SN2) reactions [START_REF] Kravtsov | Prediction of the preferable mechanism of nucleophilic substitution at saturated carbon atom and prognosis of S N 1 rate constants by means of QSPR[END_REF].

Difference descriptors

The last approach that could be used for generation of reaction transformation descriptors without explicit consideration of reaction center is based on arithmetic difference between vectors of reagents and products. The approach is based on the idea that if reagents and products are represented by fragment descriptors, difference between them represents fragments that exist in one part of the equation but absent in the other.

In such a reaction center is implicitly represented without AAM.

In case of binary fingerprint usage there are two options how the difference fingerprint could be calculated. First, the difference could be calculated as element-wise subtraction, in this case final fingerprint will contain "1", "0" and "-1" and no longer a binary number. This approach was used to classify enzymes according to reactions they catalyze [START_REF] Hu | Assignment of EC Numbers to Enzymatic Reactions with Reaction Difference Fingerprints[END_REF]. The other option is to calculate difference fingerprint as bitwise logic OR operation, resulting fingerprint will be still binary but will coincide for forward and backward reaction [START_REF]Fingerprints -Screening and Similarity[END_REF].

The drawback of the approach is that only fragment descriptors can be used (otherwise their meaning is unclear), and reaction should be perfectly balanced. If one reagent of product is absent descriptor vector will contain meaningless values. Schneider et al. [100] implemented the idea of reaction difference fingerprints, using the following procedure: authors calculated the bitstrings for each molecule involved in the reaction and the descriptor representation of the reaction was found as the difference between the sum of product bitstrings and the sum of reactant bitstrings.

Resulted vector was summed with descriptors of additives to enhance the ability of the model. It allows classifying reactions by their types.

Another type of difference descriptors has been proposed by Ridder and Wagener [101] who calculated the reaction descriptors as the difference between the frequencies of occurrence of some fragment in the products and in the reactants. A similar approach, but using fragment descriptors, which were called atomic signatures (sort of augmented atoms), was applied by Faulon et al [102]. The MOLMAP reaction descriptors calculated as the difference between the product and reactant MOLMAP descriptors, which are calculated as Kohonen's map of bond centered descriptors for a given molecule. They were used several times for classification tasks [START_REF] Zhang | Structure-Based Classification of Chemical Reactions without Assignment of Reaction Centers[END_REF][START_REF] Marcou | Expert System for Predicting Reaction Conditions: The Michael Reaction Case[END_REF][START_REF] Latino | Genome-scale classification of metabolic reactions and assignment of EC numbers with self-organizing maps[END_REF] involving chemical reactions.

CGR descriptors

Explicit representation of reaction center could overcome many problems related with descriptor generation: first, making emphasis on reaction center potentially enhance predictive ability and reduce chance to find spurious correlations; second, this method is insensitive to molecule ordering and reaction balancing. In LFER approach atoms belong to reaction center were usually selected manually choosing appropriate substituent descriptor or calculating some quantum chemical property of atoms causing interest.

Reaction center graph representation could be used as a basis for descriptor calculation. Varnek [START_REF] Varnek | Substructural fragments: an universal language to encode reactions, molecular and supramolecular structures[END_REF] noticed that CGR can be regarded as a pseudo-molecular graph and its fragment descriptors could be readily calculated using approaches designed for application to molecules, Figure 15. Resulted descriptor vectors for chemical reactions could be compared using wellknown similarity metrics and applied in reaction similarity search [103]. De Luca et al

[104] used neighborhood behavior approach [105] with different types of CGR fragment descriptors and similarity indices to find optimal strategies to assess reaction similarity.

The paper was based on the idea that the most similar reactions should correspond to the same type. Best descriptors selected by neighborhood behavior were used for mapping reactions using Kohonen's self-organizing map approach and it was shown that reactions of a given type fall into the same or nearby nodes. The latter was used to select reaction signatures corresponding to a given reaction center. Application of fragment descriptors for CGR in QSRR modeling will be described below.

Reaction condition descriptors

Reaction conditions are essential for prediction of reaction properties. Temperature, pressure, reactant concentration and solvent are the most important of them. Usually concentration independent characteristics are modeled (reaction rates, equilibrium constants, usually taken as logarithms) and thus reagents initial concentrations are neglected. Temperature and pressure are numerical values that could be easily added to the descriptor vector. Due to Arrhenius equation (Eq. 4) if logarithm of reaction rate is modelled it is more natural to introduce temperature descriptor as inverse of absolute temperature (1/T, where T is in Kelvin).

Traditionally, the mechanism of solvent action is associated with its ability to stabilize polar molecules in solution (called polarity), participation in H-bonds with solvent being H-donor or H-acceptor (called H-bond acidity or basicity) [START_REF] Reichardt | Solvents and Solvent Effects in Organic Chemistry[END_REF]. Moreover, a dissolved molecule could polarize solvent itself that changes the polarity of the latter. In modeling, solvent is represented as a vector of usual molecular descriptors (e.g., fragment descriptors and others) and by physicochemical parameters describing the main mechanism of solvent action on solute. In the work [START_REF] Kravtsov | Prediction of Rate Constants of SN2 Reactions by the Multicomponent QSPR Method[END_REF] authors included both type of descriptors however fast stepwise multiple linear regression used for descriptor selection picked only physicochemical parameters of solvents, that indirectly shows that the latter could better reflect solvent influence.

Among physicochemical parameters the most important are dielectric permittivity and refraction index. The former is related to solvent polarity and the latter with its of some molecules in different solvents.

Computational approaches to reactivity modeling

One of the most universal approaches for modeling properties of chemical objects is based on application of machine learning tools to chemical systems. Generally this approach is called QSAR (Quantitative Structure-Activity Relationship) or QSPR (Quantitative Structure-Property Relationship), with QSAR being the most widespread term. The essence of the approach is based on finding mapping between a space represented by set of structural descriptors {𝑋 ; , 𝑋 = , … 𝑋 C } and predicted property:

𝑌 = 𝑓(𝑋 ; , 𝑋 = , … 𝑋 C )
This mapping is done in a way valid for most of objects. Thus model obtained on some limited portion of data called training set could be applied for the rest of objects.

Linear Free Energy Relationships [115] use very similar technique in general but

interpretability of the model is a cornerstone of the approach. That is why in LFER linear Linear free energy approach dates back to 1937 when Hammett [START_REF] Hammett | The Effect of Structure upon the Reactions of Organic Compounds. Benzene Derivatives[END_REF] proposed it to predict acidities of chemical compounds. The name arises because the logarithm of an equilibrium constant (at constant temperature and pressure) is proportional to a standard free energy (Gibbs energy) change, and the logarithm of a rate constant is a linear function of the free energy (Gibbs energy) of activation. The approach was extensively used to establish linear correlations between substituent or solvent descriptors and reaction characteristics such as reaction rate or equilibrium constants. There were a lot of investigations done in this field that are reviewed in the book of Palm [116]. In the end of XX century with the development of computers quantum chemical approaches have mostly substituted LFERs. Usually quantum chemical studies involve building simple correlations between parameters of molecules and their reactivities which are alike LFER, for example, linear equations between Diels-Alder reaction rate and conceptual DFT indices [117]. Reaction kinetics study sometimes still involves LFER to understand substituent or solvent effect and thus get some ideas on reaction mechanism.

In this part we will describe the application of different approaches to study chemical reactions.

Quantum chemical calculations

The quantum chemical methods are based on approximate solution of Schrodinger Energy and its derivatives are used to locate extremums on potential energy surface corresponding to stable chemical structures or transition states and in such a way used to assess molecular energy and geometry. Many other characteristics of molecular systems (including spectra) could be calculated in a similar manner.

Predicting of thermodynamic parameters

The quantum chemical calculations are the main way to predict thermodynamic parameters of molecules nowadays. Usage of quantum chemical approaches for prediction of different thermodynamic parameters is published each year in thousands of articles. Special approaches for precise estimation of thermodynamic properties of molecules taking into account all required energetic terms and usually some empirical corrections are proposed [127,128].

Prediction of thermodynamics parameters for reagents and products could be used to assess reaction enthalpy and free energy. The latter could be used for prediction of reaction equilibrium constant. For example, the company "Schrodinger" proposed an accurate quantum chemical approach for assessing tautomeric equilibrium constants [129].

Despite quantum chemical approaches are the most used methods for prediction of reaction thermodynamics, they are rather slow and the need to estimate solvent effect sufficiently complicates this problem. The achievement of chemical accuracy in the prediction requires the performance of very sophisticated and resource-intensive calculations. And even with them the quality of data is far away from experimental accuracy [START_REF] Peverati | Quest for a universal density functional: the accuracy of density functionals across a broad spectrum of databases in chemistry and physics[END_REF]. Results of a blind competition performed in the frame of SAMPL2

Challenge [130] where various, mainly quantum chemical, approaches competed in the prediction of solvation energies and energy differences of tautomers have shown that computational accuracy they reach now (the mean square error of the energy calculation was ~2.5-3 kcal/mol) is insufficient for the adequate description of the energy and is much lower than the experimental one. Almost all methods used had fitting parameters adjusted to a training set proposed in the competition that allows more appropriate estimation of desired property compared to ab initio calculation.

Mechanism of reaction

Usually mechanism of reaction is explained in terms of Transition State Theory, that was developed simultaneously by Eyring [131], Evans and Polanyi in 1935 [132]. A minimal energy pathway on potential energy surface that connects reagents and products is stated in the theory. This pathway describes alterations of geometry and energy in the system while reagents are transformed to products. Minimal energy pathway is the pathway that requires the least amount of energy. The movement along the pathway requires changing of one or more molecular coordinate and is called reaction coordinate. Eyring have shown that reaction rate depends on activation free energy:

Eq. 5. 𝑘 = ^_à 𝑒𝑥𝑝 c- ∆d e ?`f
where ∆𝐺 g is known as the Gibbs free energy of activation, is the standard molar Gibbs energy change for the conversion of reactants into transition state. A plot of standard molar Gibbs energy against a reaction coordinate is known as a Gibbs-energy profile. In principle the right-hand side of equation (Eq. 5) should be multiplied by a transmission coefficient, κ, which represents the probability that an activated complex forms a particular set of products rather than reverting to reactants or forming alternative products.

In such a way, quantum chemical calculations could be used to understand mechanism of reaction and to find its rate constant [133]. Even for simple reactions quantum-chemical study of its mechanism and kinetics is computationally intensive, but nevertheless it is the most developed and widely used approach for reaction study and qualitative assessment of their characteristics [134].

QSAR in reaction modeling

Chemoinformatics approaches were rarely applied to study and predict reaction characteristics. Mainly it is related with reaction complexity and lack of data ready for modeling. Thus works related to reaction modeling appeared sporadically and no systematic studies were still published.

Reaction kinetics

The collection of data on reaction kinetics were mainly inspired by great interest in reaction mechanism study and popularity of linear free energy relationships in the XX century [116,135,136]. By the end of 1990 th general workflow of QSAR modeling had been developed and thus first attempts to build models for reaction kinetics based on QSAR methodology was done.

The first works in reaction rate prediction were done in N. Zefirov's group and were based on application of artificial neural networks and molecular descriptors of different type to chemical reactions. In the work of Sukhachev et al [START_REF] Sukhachev | QSPR approach to the calculation of rate constants of homolysis of nitro compounds in different states of aggregation[END_REF] rate constants of nitro compounds decomposition in gas phase were predicted. Authors used ordinary topological and fragment descriptors ignoring conditions that were the same for all measurements.

In the work of Halberstam et al [137] a model for the rate constant of 2092 acid ester hydrolysis reactions was built. Since one of reagent was unchanged (water) the descriptors of chemical transformation were simply molecular descriptors of ester, quantum-chemical descriptors were used. To represent reaction conditions reaction temperature and Palms descriptors of solvents were concatenated with structural descriptors. On validation set (209 reactions), RMSE of prediction was 0.34 logk units.

Using the same dataset and workflow but fragment descriptors to represent ester structure Zhokhova at al [138] slightly increased the quality of prediction-RMSE dropped down to 0.31 logk units.

For bimolecular nucleophilic substitution reaction rate assessment Kravtsov et al [START_REF] Kravtsov | Prediction of Rate Constants of SN2 Reactions by the Multicomponent QSPR Method[END_REF] proposed approach where molecular descriptors of two substrate and product molecules were combined into the feature vector using for modeling by means of neural

Reaction conditions prediction

Prediction of optimal conditions for chemical reactions is very interesting but yet really insufficiently explored topic in reaction modeling. One of the first work in this direction was the one by Struebing et al [START_REF] Struebing | Computer-aided molecular design of solvents for accelerated reaction kinetics[END_REF] where optimal solvents for nucleophilic substitution reaction were predicted. Proposed approach used quantum chemical calculations of transition states of reactions in a given solvent with the following LFERlike equation that united activation energy and solvent parameters. The latter was used to find solvent that potentially enhances reaction rate and the procedure repeated until selfconsistency reached. As the result authors reported about new found solvent for reaction that increased reaction rate on 40%. This approach was quite resource-consuming -one prediction required some days of calculation.

Chemoinformatics approach in this case could be much faster. Marcou et al [START_REF] Marcou | Expert System for Predicting Reaction Conditions: The Michael Reaction Case[END_REF] proposed a classification model that could predict general type of catalyst and solvent for Michael reactions. Recently, approach to assess optimal conditions for hydrogenation reactions was proposed by Lin et al [149]. The training set reactions were extracted from Reaxys database and curated in fully automatic manner. Optimal conditions were predicted using similarity based approach: for test set reaction system looks for most similar reaction center environment in curated database of reactions where required transformation proceeds and where does not. The approach was validated internally and externally and has shown very good results.

Reaction yield prediction

Prediction of expected yield for reaction running under certain condition could be rather interesting for synthetic chemist to assess yield of multistep procedure or ensure that reaction conditions were selected properly. There are two very recent publications that dealt with yield prediction.

Very recently Skoraczynśki et al [150] modeled 425 000 of reactions from Reaxys. In this chapter technologies used in the work for reaction modeling will be reviewed.

Moreover, some clarification of technical aspects of the workflows used in the work is given.

3.1.

Quantitative Structure Reactivity Relationship (QSRR) modeling strategy QSAR/QSPR methodology application to chemical reaction proposed and used in the work will be described in the following chapters. In order to discriminate it from other approaches we called it for convenience Quantitative Structure Reactivity Relationship (QSRR).

The general workflow of QSRR is the same as for QSAR study and starts with generation of descriptors for the objects of the dataset. The next step is selection of model validation procedure and machine learning method. Hyperparameters of machine learning method and fragment descriptors type maximizing model performance are selected using stochastic algorithms. The parameters shown best performance in cross validation are used for final consensus model building. All the steps will be discussed in chapters below.

Reactions Descriptors

Every modeling procedure starts from descriptors generation since machine learning methods usually require input represented as numerical vector. In this work two types of fragment descriptors were used for modeling of reactions -SiRMS and ISIDA.

The SiRMS descriptors were used for modeling of reaction for the first time in this work. The procedure of generation descriptors is fully described in the article on E2 reactions modeling in the chapter 6.2, as the procedure that was used only in that part of the work and not included in general workflow. SiRMS descriptors were calculated using • Minimal length from 2 to 4;

• Maximal length from 3 to 8;

• Fragment types: sequences and augmented atom type descriptors including description of only atoms, only bonds or both atom and bonds;

• Formal charge on atoms were optionally explicitly included into fragment description;

• Shortest paths or all possible paths were optionally explored in sequence descriptor generation.

As a result, 616 fragmentation schemes were used. refractive index is denoted as n for simplicity). In order to describe water-organic solvent mixtures descriptor specifying the molar ratio of organic solvent in water was added. The inverse absolute temperature, 1/T (in Kelvin degrees) was also used as temperature descriptor.

Model Validation

The best parameters for machine learning methods are selected within 5-fold cross-

Machine-Learning methods

Machine learning algorithms can fit some flexible function to describe the data, and thus it could extract rules from observed data without being explicitly programmed [START_REF] Varnek | Chemoinformatics as a Theoretical Chemistry Discipline[END_REF].

Such algorithms overcome following strictly static program instructions by making datadriven predictions or decisions, based on sample inputs. Machine learning is widely used in a range of computing tasks where designing and programming explicit algorithms with good performance is difficult or infeasible.

Support Vector Machine

Among the most popular state-of-the-art algorithms is Support Vector Machine (SVM) -a supervised non-probabilistic learning method used for classification (called SVC or simply SVM) and regression analysis (SVR).

In classification task a set of training examples is given, each labelled as belonging to one of two classes, an SVM training algorithm assigns new objects some class. An SVM looks for hypersurface in some descriptor space that separate objects of different classes.

The surface is adjusted in a way that the gap (called margin) between objects of opposite classes is as wide as possible. New examples are then mapped into that same space and predicted to belong to a category based on definite side of the surface they fall on.

In addition to performing linear classification, SVMs can efficiently perform a nonlinear classification using so called kernel trick, implicitly mapping their inputs into highdimensional feature spaces.

Consider some training data 𝐷 that are represented by a set of n objects with known class:

𝐷 = {(𝒙 𝒊 , 𝑦 j ) | 𝒙 𝒊 ∈ ℝ n , 𝑦 j ∈ {-1, 1}} jp; C
where the 𝑦 j is either 1 or -1, indicating the class to which the point 𝒙 𝒊 belongs.

Each 𝒙 𝒊 is a pdimensional real vector. Optimal classifier is found by selection of the maximum-margin hyperplane that divides the points having 𝑦 j = 1 from those having Thus the following task of quadratic programming is posed: where C is the cost that define the severity of the penalty for points whose predicted value deviates from experimental one more than ε, Figure 18: Kernel trick based implemented using dual formulation of SVR optimization task is utilized to build non-linear regression model.

𝑦 j = -1,
; = 𝒘 `𝒘 → min 𝑦 j (𝒘 • 𝑥 j -𝑏) ≥ 1
𝑦 j -(𝒘 • 𝒙 𝒊 ) -𝑏 ≤ 𝜀 + 𝜉 j (𝒘 • 𝒙 𝒊 ) + 𝑏 -𝑦 j ≤ 𝜀 + 𝜉 j * 𝜉 j ≥ 0 𝜉 j * ≥ 0

Random Forests

Another popular machine learning method is Random Forest (RF). Random forests are a combination of tree predictors such that each tree is built on random subset of objects sampled independently and with the same distribution for all trees in the forest [153]. This method is evolution of simple decision tree [154], where each node of tree is a rule based on value of some descriptor that selects subset enriched with objects of given class. Thus a tree is an ensemble of rules that leads to decision about class attribution of an object.

In random forest tree predictors are united in a way that every tree gives independent prediction and final decision about class attribution is based on simple majority of votes. For random forest regressor trees for prediction of continuous variables are used and averaging predicted values of trees makes decision. Each tree is built on subset of objects selected by bootstrapping [155]. Moreover, during tree learning (growing) procedure, when algorithm looks for feature giving best split in branching, random subset of descriptors is considered rather than all descriptors. The main hyperparameters in RF is number of trees and number of descriptors in random subset.

The more trees are generally better for this machine learning method, but it also leads to increase in resources and time consummation. Optimal number of descriptors in random subset is usually adjusted using cross-validation or on the basis out-of-bag sample prediction.

GTM

As a visualization tool Generative Topographic Mapping (GTM) was used. In GTM, each point in the low dimensional (usually 2D) latent space (LS) is mapped onto the manifold embedded in a high-dimensional data space (input space, IS), Figure 19.

The manifold is defined by a mapping function y (x; W) assessed with the help of m radial basis functions (RBFs) of width w regularly distributed in LS. The latent space is covered by a mesh containing k nodes each of which corresponding to a normal probability distribution (NPD) centered on the manifold in IS. Ensemble of NPDs is used to compute a posterior probability for a data point tn in D-dimensional IS to be projected onto a node xk:

𝑝(𝑥 ^|𝑡 C , 𝐖, σ) = p(𝑡 C |𝑥 ^, 𝐖, σ) ∑ p ^(𝑡 C |𝑥 ^, 𝐖, σ)
where W is a parameter matrix and σ the variance of the distribution of t:

𝑝(𝐭|𝐱) = ‰ 1 2𝜋𝜎 OE • = exp •- 1 2𝜎 ‖𝐲(𝐱; 𝐖) -𝐭‖ = "
The log likelihood of the whole data set is calculated according to equitation:

ℒ(𝐖, 𝜎) = | 𝑙𝑛 - 1 𝐾 | 𝑝(𝑡 C |𝑥 ^, 𝐖, 𝜎) ^- C
The GTM is optimized with an expectation-maximization (EM) algorithm using data likelihood (ℒ) as the objective function (the best GTM map corresponds to the highest ℒ). The mapping depends on four parameters: the number m of RBFs, the grid resolution k, the RBF width w, and the weight regularization coefficient l. The latter is used for re-estimating the W parameter matrix and influences the flexibility of the manifold. Notice that a too flexible manifold, although nicely approximating the training data, may lead to overfitting when training set objects ideally fit manifold but new data points are located too far from it. In this work GTM map parameters w and l were optimized by likelihood.

GTM approach is unsupervised machine learning methods, as it does not use target values while learning. Obtained map however can be colored by different properties to define zones of selectivity on the map, as it was done in our previous work [156]. This approach is used in this work for the analysis of SN2 reactions chemical space.

Applicability domain

All QSAR models are strongly connected with the training set. Due to limitated number of trained instances prediction could be done more or less confidently for objects similar to training set ones. The set of objects that could be predicted confidently by a model is called its applicability domain (AD) [157], Figure 20. The test set compound 1 (in green) is inside the AD and its prediction is considered reliable while the test set compound 2 (in red) is outside and therefore, its prediction is considered unreliable. The concept of applicability domain is extensively studied in chemoinformatics [157][158][159][160][161] since datasets are usually rather limited and machine learning methods being good in interpolation are unstable in extrapolation.

Genetic Algorithm driven parameters optimization

The fragment descriptor types (fragmentation scheme) and hyperparameters for machine-learning methods were selected within stochastic search approach. As an algorithm was selected program by Horvath [162] with genetic like algorithm of search.

Predictive performance of models based on given vectors of parameters to be optimized is checked within cross-validation procedure. The first generation vectors are initiated unsystematically. Two sources of new vectors come from this time : randomly generated and generated by crossover using best-performing vectors. This strategy helps to avoid full iteration of all possible variants and gives reasonably good results in acceptable time.

Model publishing

Modeling procedure was made by in-house program called ChemoInformatics and Molecular Modeling Lab tools (CIMMtools). This program produce special compressed container with model. Such containers are stored at our server and are used by special dispatcher. All obtained models are available at http://cimm.kpfu.ru/predictions -visual interface for working with model dispatcher. Moreover, server have special module, that takes input reactions and makes all preprocessing before prediction, according to specification of model that will be used. Then model gives prediction for curated reaction and gives evaluation of prediction by bounding box applicability domain. The server can handle several users simultaneously or make a queue of queries. Usage of models is free of charge, but requires registration procedure on the site. The example of model application is given in chapter 4.4.

Matched Molecular Pairs

Matched Molecular Pairs (MMP) have generally been defined as "molecules that differ only by a particular, well-defined, structural transformation" [163]. The analysis of such pair is usually used by cheminformatics to see smoothness in change of the properties of two molecules that differ only by a single chemical transformation, such as the substitution of a hydrogen atom by a chlorine one (Figure 21). Since the structural difference between the two molecules is small, any experimentally observed change in a physical or biological property between the matched molecular pair is attributed to substituent effect.

Quantum Chemistry calculations

Density functional theory (DFT) is a computational quantum mechanical modelling method used in physics, chemistry and materials science to investigate the electronic structure (principally the ground state) of many-body systems, in particular atoms, molecules, and the condensed phases [164]. DFT states that ground state energy of the system could be expressed as the functional of the electron density unlike classic quantum chemistry (Hartree-Fock and post-Hartree-Fock methods), which relates energy with the many-body wavefunction. Whereas the many-body electronic wavefunction is a function of 3N variables (the coordinates of all N particles in the system) the electron density is only a function of x, y, z -only three variables. To calculate energy of the molecular system and other properties one need to know exchange-correlation functional that approximates dependence of exchange and correlation energy on electronic density or so called Kohn-Sham orbitals.

In this work for DFT calculations Priroda [165] program was used, that is one of the fastest DFT code due to approximation of Coulomb and exchange-correlation terms After validation procedure, new window with AAM established by embedded algorithms is given to user to check it and modify if needed. On the same page one can select the model to apply and choose conditions (Figure 24). If single molecule is entered, reaction models are not available.

Figure 24. Window with selection of model and conditions for modeling.

After the selection of all parameters and clicking "Modeling" button, the server returns the result of modeling (Figure 25). At this window predicted value and variance of prediction within consensus is given. Also results of AD application are given in line "Trust of prediction". The trustworthiness of prediction depends on fragment control applicability domain and variance of prediction. If prediction of test set appears to be out of fragment control AD for more than certain percent of models inside consensus or variance exceed certain value, trustworthiness will be lowered from "optimal" to "good".

If two conditions are fulfilled simultaneously it becomes "poor". The reason of "Trust of prediction" lowering is given below in "Distrust reason" section.

Figure 25. Results of modeling procedure.

Hense, modeling procedure is user-friendly and includes 3 simple steps: inputvalidation -results. All the preprocessing and descriptors generation is written in the model itself and hidden from the user. As a result fast and easy tool for reaction modeling was obtained.

In-house software tools

All the tools developed in the work are written in Python 3.

Development of a comprehensive reaction database

Nowadays, modern chemical databases (Reaxys, CASReact, InfoChem, etc) contain some 100 million reactions. They store every reaction as a separate record that consist of reagents, products and some additional information. Despite the fact that amount of known reactions is incredibly high, there are almost no data for modeling reaction rates, almost no data on other kinetics parameters annotated in databases. Usually only reaction time and yield of desired product are annotated according to the needs of synthetic chemists. Recent work shows that these parameters could hardly be modelled [150]. Thus the content of existing reaction databases is not satisfactory for reaction characteristics prediction.

From the other hand, standardization of reaction is tricky and not so well elaborated as for molecules. The reactions are usually represented in the way they were in original article, databases are just storage of information about successful reactions from the article with no moderation. Some preprocessing is done to merge reactions with the same structures of reagents-products. Notice that reactions with the same transformation and conditions are usually saved. Usually reactions in such databases are stochiometrically unbalanced (some reactants and products are lost), its additives and conditions are represented in text fields in non-standardized way and could be lost, for example, only half of reaction extracted from Reaxys had temperature in corresponding field [149].

To overcome the problem of lack of data, a new database of kinetic measurements of reactions was collected. Relational database was designed to deal with big amount of data on different reactions that can differ by transformation or/and conditions.

In this chapter the workflow of standardization, storage, representation and filtering of reactions will be described. This workflow appeared as the result of further development of CGR concept and its representation.

CGR technology development

The new challenges were come across during the realization of new database and standardization techniques that are related to the needs in further development of CGR technology, its adaptation to the challenges we faced upon reaction modeling and storage.

In this chapter the developments in CGR technology such as introduction of dynamic atom concept, more useful way of CGR storage and creation of reaction center signatures based on CGR are described. Such features are needed to increase the speed and easiness in processing and storage of reactions.

CGR technology extension

In this chapter the reasons for dynamic atom concept introduction into CGR approach is described.

Usually for representation of molecules and reactions reduced graph representation with implicit hydrogens is used. It assists not only to save data storage space, but reduce number of fragments descriptors that will be generated for the structure. Usually it does not cause problems, since the change of hydrogen position is encoded in change of bonds orders of between heavy atoms. But in some cases, e.g. for zwitter-ion formation reactions (type of tautomeric equilibria) or SN2 reactions with OH-or NH-containing nucleophiles, implicit hydrogen migration is not captured in CGR since bonds between heavy atoms are not touched. Thus, only charge of some heavy atom is changed in reaction. There is several ways to encode hydrogen migration in CGR:

• Explicit hydrogen representation. Dynamic bond will appear between heavy atom and proton and thus migration of hydrogen will be encoded, Figure 26. From the other side, it will increase number of descriptors and the time required for its calculation.

• Usage of pseudoatoms representing atomic charge. The idea was to represent additional features of atom as connections from it to pseudo atoms standing for charge label (+1/-1/0), stereo label (R/S), etc. It was proposed by Jauffret [START_REF] Jauffret | Machine learning of generic reactions: 2. toward an advanced computer representation of chemical reactions[END_REF]. Migration of hydrogen will be encoded as dynamic bonds between heavy atom and virtual atom that represent charge, Figure 26. However, it causes a problem with visualization and CGR could hardly be represented by widespread SDF format. Additional atoms will also drastically increase number of descriptors that decrease speed of model building.

• Dynamic atom. The third idea was to introduce analog of dynamic bond for one sole atom -dynamic atom. Dynamic atom encodes changes in some atomic property (charge, stereochemistry) in chemical reaction, Figure 26. This will solve a problem with encoding proton migration, as atomic formal charge is changed that could be captured by dynamic atom label. Generally the number of descriptors does not increase, as dynamic charge is a simple atomic label. This label can be correctly represented in SDF format and visualized by common visualization tools (e.g. ChemAxon Marvin View). Thus, this approach was implemented in our SDF specification.

CGR storage in SDF format

In this chapter reasons of SDF format usage for CGR storage will be discussed and features allowing to do it natively for SDF file format will be described.

MDL/SDF [START_REF] Dalby | Description of several chemical structure file formats used by computer programs developed at Molecular Design Limited[END_REF] is one of the most popular formats in the world for storage information about molecules. (starting with "M SDT") contains title of property: "dynatom" and "dynbond" for dynamic atom and bond correspondingly. The third line starting with "M SDD" contains information needed for visualization of CGR -coordinates on the screen where to place the label. Forth line starting with "M SED" contains value of dynamic atom and bond label. Dynamic bond labels representing changes in bond orders are given using the following scheme "R>P", where R and P bond orders of the bond in reagents and products (0 for absence of bond). Dynamic charge is represented using "c+1" and "c-1" labels standing for charge increase or decrease by 1 correspondingly. The dynamic labels were proved to be rendered by ChemAxon programs and the latter could be used for creating CGRs in proposed specification.

CGR signatures

At this chapter such newly developed features as canonical CGR and reaction center signatures will be described. These two techniques is used for fast duplicated transformation search and reaction classification without standard graph embedding procedure.

For canonical numbering Morgan [182] like algorithm was developed based on prime number usage as in paper of Ihlenfeldt and Gasteiger [183]. In it we use the property to that every number could be factorized uniquely. Unlike classic Morgan algorithm where extended connectivity index for a given atom is calculated summing up corresponding indices of surrounding atoms, we use multiplication of prime numbers representing atomic connectivity. In this case, there is no possibility that two atoms have the same extended connectivity index by chance.

Conceptually algorithm includes several steps. The first step includes ranking of tuples describing atomic properties (element, connectivity, charge, isotope, stereolabel, etc) of every atom in molecule as in CANGEN algorithm for SMILES canonization [START_REF] Weininger | SMILES. 2. Algorithm for generation of unique SMILES notation[END_REF].

Based on the rank obtained prime number is selected from prime number table. Thus every atom assigned prime number that coincides for atoms having the same type. This number used to define amount of unique atoms. At the next step this number for a given atom is multiplied to the corresponding numbers of its neighbors and again used to define amount of unique atoms. Obtained numbers are sorted and again used for selecting prime number from table. This procedure is repeated several times, until the number of unique atoms will not change for 3 following cycles. This final numbers are sorted and the ranks are used as canonic atom numbers. The latter are applied for linear string generation using rules similar to canonic SMILES generation [START_REF] Weininger | SMILES. 2. Algorithm for generation of unique SMILES notation[END_REF].

For CGR generation current implementation generates canonic numbers for all molecules in reagent and product side and SMIRKS-like string is generated using All described signatures are quite long strings thus we compressed them using hash function and represented as hexadecimal number. Developed signatures help to match reactions superfast without expensive graph embedding procedure, cluster them and analyze influence of reaction center surroundings.

Data collection

As we already mentioned there exists no database of well-structured data on chemical reactions, especially its kinetic characteristics that are ready for modeling. We decided to collect our own dataset of kinetic and thermodynamic properties of chemical reactions. Data were collected manually from the reference book by Palm [START_REF] Palm | Tables of Rate and Equilibrium Constants of Heterolytic Organic Reactions[END_REF] and PhD theses that were defended in Kazan Federal University, Russia. The information was extracted manually. In this chapter the rules that were used for selection and annotation of data will be described. These rules were used to minimize efforts needed for data standardization; some of them were introduced to make possible usage of the other modeling approaches rather than CGR used in this study.

Datasets of SN2 and E2 reaction rate and equilibrium constants of tautomerization was collected from book edited by Palm [START_REF] Palm | Tables of Rate and Equilibrium Constants of Heterolytic Organic Reactions[END_REF], cycloaddition (CA) reactions rate constants, activation energy and pre-exponential factor were extracted from PhD theses of Kazan Federal University (totally about 20 theses). The databases were created using InstantJChem [187] database management system.

We posed some constraints on reaction to be extracted. For SN2, E2 and CA reactions, only constants that were declared for bimolecular processes were added to database. Only reactions in pure solvent or water-organic solvent with known molar ratio of solvents were annotated. In all cases reactions with additives and catalysts were ignored, as well as reactions under not standard pressure (other than 1 atm).

The rules of reaction annotation were slightly different for different dataset. Order of molecules in reaction was fixed. Molecules in SN2 reaction were drawn in following order: substrate, nucleophile, anion, cation -> major product, minor product, anion, cation. For E2 the order was following: substrate, base/reagent, counter ion of base/reagent -> product, cleaved group, changed base/reagent, counter ion of base/reagent (see Scheme 2). For cycloaddition reaction diene was drawn first, dienophile second. Counter cations or anions in all type of reactions are drawn after reagents and products.

Scheme 2.

Usually description of reaction in literature sources included only description of reagents. It is not problematic for most of reaction types to guess the product, however for CA reaction it caused unclear stereoconfiguration or even regioisomer of products, Figure 28. Consequently, CA reaction data was stored as a mixture of all possible reactions with special markers to identify them while dataset preparation. The tautomeric equilibria database has its own peculiarities. Only equilibria between two possible tautomers were collected. Since equilibrium is reversible reaction it is important to select which molecules is reagent and which is product. The roles of molecules were assigned according to the type of tautomerism, e.g. for keto-enol equilibria keton was reagent and enol was product.

For all databases ionic compounds were represented as two separate ions. Bonds of metals with heteroatoms -N, P, O, S, Se, Te, F, Cl, Br, I atoms -were considered ionic except of organometallic compounds. All molecules were stored as graphs with implicit hydrogen. All known information about reaction was extracted.

Mentioned rules helped to collect database that have unified representation of reaction. The amount of collected data for each reaction type is given in Chapter 4.1.5.

Data curation

Initial data set should be checked for possible errors in representation or in values of the property. This chapter will describe techniques of data curation and stages of this process. This workflow was used for preparation the dataset for modeling.

Quality of data is a very important issue for building predictive models using machine learning tools. While best practices for data curation in QSAR modeling [188,189] and chemogenomics data cleaning [190,191] were published there is no such kind of information for chemical reactions. Reaction data curation is more challenging task than molecular data curation mainly due to the fact that conditions should be taken into account. From the other side Arrhenius law could be used for finding doubtful data on reaction rate. Thus, we developed our own workflow for processing chemical reaction that takes into account the features of reactions as chemical objects.

Reaction cleaning process includes two stages: ( 1 Usually the same transformation could proceed under different conditions and thus with different rate constant. For further usage in modeling it is more useful to group all reactions having the same transformation. Then dataset could be stored as relational table.

Structure standardization

The first step is standardization of reaction representation is unification. In order to find group of reactions with same transformation that differ in conditions, initial database structures was standardized in order to remove different representations of the same molecules using ChemAxon Standardizer [192]. were not mapped by mistake).

Atom-to-atom mapping

Second step, that will be described in this chapter, is devoted to search and correction of AAM by means of new CGR feature as linear representation of CGR and canonical SMIRKS.

CGRs could be used to identify incorrect AAM. CGR corresponding to incorrect AAM usually have more dynamic atoms and bonds than the correct one and thus is more complicated that is consistent with minimum chemical distance principle [193], Figure 29.

For every transformation we created two alphanumerical strings: canonical SMIRKS representation of reaction ignoring atom-to-atom mapping (CSMIRKS) and CGR signature (CGRS). The latter represents a unique identifier of reaction transformation like SMILES or InChI for molecules. Both line notation strings were created by in-house program that uses aforementioned Morgan-like algorithm for canonical numbering. Then canonical name is created according to defined numbering using SMILES-like rules. For human readability CSMIRKS and CGRS long text strings of variable length are compressed using MD5 hash function [194] and represented as hexadecimal number that are called CSMIRKS-key and CGRS-key correspondingly.

Transformation with the same set of reagent and products will have same CSMIRKS whichever AAM is. CGRS will depend on AAM: two formally different but correct AAM (reactions A and B on Figure 29) have the same CGRS while they do not coincide for reactions with correct and wrong AAM (reaction A and C on Figure 29). As there were no reliable and numerous data on rates of reaction involving stereoisomeric reagents or with formation of isomeric products it was decided to ignore stereochemistry of reaction and to use 2D descriptors. Thus upon CSMIRKS and CGRS generation stereochemistry of reaction was ignored and transformations with different absolute stereolabels of reagents/products were considered the same. Reactions having the same CSMIRKS are grouped together. Different CGRS for reactions having the same CSMIRKS could be caused only by AAM error. Such reactions were identified and their AAMs were manually fixed.

The described approach identifies AAM errors for duplicated transformations however it does not work if the transformation is performed in one sole condition or the error was made for all transformation of given type. Thus, we used approach for checking validity of AAM by verifying signatures of reaction centers.

Reaction center is a set of atoms that change their environment in reaction.

Generally, complex reactions could have more than one reaction center and thus deletion of conventional atoms and bonds will lead to disconnected graph. Thus, reaction center is a one connected component of such graph. Reaction center signatures for a given transformations were generated using approach described earlier in chapter 4.1.1.3. For a given reaction type the number of signatures is limited, e.g. for SN2 reactions cleaned dataset it is equal to only 29. Error in AAM lead to appearance of wrong reaction center (Figure 29). They could be easily manually distinguished and AAM of corresponding reactions were fixed.

Curation of temperature data

The next task is curation of reaction conditions. The most important conditions for us is solvent and temperature. Error in solvent specification could be corrected only by human and could hardly be identified automatically. Temperature data being the most important characteristics could be automatically verified.

Since we mainly interested in reactions between dissolved molecules, solvent or mixture have to be in liquid phase. All datapoints with temperatures below freezing point or higher than boiling point of pure solvent or component in the solvent mixture were discarded as doubtful.

As a second filtering rule we required that temperature dependence should follow chemically meaningful behavior. Since the activation energy (Ea) is unknown for reaction of dataset Arrhenius equation could not be applied and Van't Hoff rule was utilized. The latter states that reaction rates usually increase 2-4 times every 10 degrees Celsius [195].

This law can be used as a filter criterion for anomaly change of rate constants with change of temperature.

As the second check, we applied Arrhenius equation to compare temperature and reaction rate constants for different reaction conditions description corresponding to the same transformation. We mark condition description of reactions as suspicious if two same transformations differ only in temperature, the difference in temperature is less than 10°С but the difference in logk is greater than 1 (according to Arrhenius equation it should be not more than 0.6). Suspicious conditions were manually examined. If there were examples supporting that the value is correct (similar measurements from other references), this datum was considered as valid. If there were not such examples, datum was discarded from modeling set.

Duplicate detection

This chapter is devoted to decisions that should be taken in final dataset preparation and how duplicated data are removed.

In this work measurements of rate constant corresponding to the same transformation considered as duplicates if they have the same conditions that are considered valuable for reaction rate constant: temperature, solvent and concentration of organic solvent in mixture. Thus, measurements taken from different literature source and records differed only by concentration of reactants were considered as duplicated reactions. The latter was done, because the second-order rate constant that was collected should not depend on concentration of components (there are some exceptional cases that are discussed below). Since reaction stereochemistry were neglected, reactions involving diastereoisomers and enantiomeric molecules running in same conditions were considered as duplicates as well.

Descriptions of conditions of duplicated reactions were combined and its rate constant logarithm (logk) were averaged if the difference in the property is lesser than 0.5 log units. Otherwise, they are marked as suspicious. The threshold 0.5 was taken from comparison of rate constant logarithm for duplicated reactions. The difference in logk for SN2 reaction set (the biggest collected set) was mainly within 0.5 log units (Figure 30). We consider this value as estimation of reproducibility of reaction rate value. Suspicious reactions were curated manually. If there were some evidences that the logk value is correct (for example, close values reaction constant of other duplicates) the value is accepted and averaged with other duplicates if any. Otherwise, values were discarded from the final modeling set. From the final dataset we excluded reactions that were performed in solvents, for which solvent descriptors described below have not been known.

The described workflow for data curation results in the curated dataset that was used for modeling.

Relational tables

Preprocessed data need to be stored in a database specifically designed for the reaction storage. This database will be briefly described in this section.

Relational database was found to be optimal for modeling dataset storage. It is logical, because database contain a lot of supporting tables about person that added record, time, molecules that contained in database, condition, properties, etc. Reaction with the same structure of reagents and products (we call it transformation) could be performed in different conditions. Thus, one-to-many relations between structure and conditions is straightforward solution for data storage which is useful for modeling. The main two tables used for reaction storage are, Figure 31: All approaches developing for data curation procedure become essential part of the database and showed their effectiveness in storage of all data that were collected manually and have some errors. The more database is flooded with data the more efficient, useful and the same time robust automatic procedures become. Each new entity that goes into database is checked by CGRS and CSMIRKS. If there is no such CGRS in database, but CSMIRKS already exist, new reaction with high probability has mapping error that needs to be fixed. Thus, we created a tool that incorporates our experience in reaction data processing and allows fast and feasible reaction search and data curation.

•

Content of the database

Summary of results of the data collection and curation is given in Table 5. Data annotation was performed by several people and took about 4 years.

In total, about 10998 records of raw data of bimolecular nucleophilic substitution (SN2), bimolecular elimination (E2), Diels-Alder (DA) reactions and tautomeric equilibria (TAU) were collected (Table 5). "Data collected" column in Table 5represents total number of reactions collected. "Curated data" shows the number of reactions after structural curation, AAM correction and deletion of full duplicates, i.e. entries where all fields coincide. The latter reason, being the major one, is caused by the fact that the dataset was collected simultaneously by several people and sometimes by error the same reaction was extracted twice. Errors in structure or absence of mandatory field were the second reasons for the deletion of datapoint. "Model set" column contains reaction selected for modeling. Data curation procedure in this case includes condition curation, duplicate identification (as described in chapter 4. 1.3.4) and averaging properties of duplicates. Total number of transformation types (combinations of reagents and products) is given in field "Transformations in model set". Difference between model set size and number of transformations reflects the fact that the property for some transformation were measured in several conditions. o molar percent of organic solvent in mixture with water (100% for pure solvent). Except of DA reactions providing in pure solvents.

• Reference -the link to the source of information:

o table in reference book [START_REF] Palm | Tables of Rate and Equilibrium Constants of Heterolytic Organic Reactions[END_REF] or thesis from which data were taken from, o page in reference book [START_REF] Palm | Tables of Rate and Equilibrium Constants of Heterolytic Organic Reactions[END_REF] or thesis.

• Property:

o Reaction rate or equilibria constant.

The data were collected using InstantJChem tool [187]. Atom-to-atom mapping was done manually according to the mechanism or reaction reported at the reference book [START_REF] Palm | Tables of Rate and Equilibrium Constants of Heterolytic Organic Reactions[END_REF]. Sets of reactions kinetics data and tautomeric equilibria were stored separately, since workflows for standardizations are different (molecules in tautomeric equilibria dataset should not be tautomerized during standardization).

The initial dataset was cleaned and converted into our own relational database format. As the result 10 903 curated reactions are stored for CA, SN2, E2 reaction and tautomeric equilibria. The developed tools and database was shown to be extremely useful for curation and storing the data.

Modeling procedure

All QSRR models that are mentioned in this chapter have the same modeling procedure. The workflow is shown on the Figure 32.

Descriptors preparation

Validation procedure and search of parameters bonds were generated. Additional options like formal charge specification and all possible path explorations were optionally used as well. In total 616 descriptors sets were generated for each dataset. Fragment descriptors vector was concatenated with descriptors of solvent, temperature, percentage of solvent in mixture with water. All resulted descriptors sets were used for model building.

• Selection of best model parameters Evaluation of each model was made within standard 5 fold cross-validation procedure that was repeated 10 times. Repetitions were made to exclude influence of random fluctuations of dataset composition on estimation of model performance. It leads to increased time for model preparation that is why checking all sets of fragment descriptors with all possible hyper parameters for SVM is impractically long. To avoid this time consuming calculation, genetic algorithm [162] was used for the best model selection. 10 best models returned after certain number of epochs that differ in kernel or/and descriptor set were selected for consensus model. The training set datasets in cross validation were used in consensus. Thus, consensus includes 500 individual models: 10 best combinations of fragment descriptor type with SVM kernel times 50 training sets used in cross validation (5 folds times 10 repetitions).

• Model performance evaluation

Test set predictions of individual models in consensus within cross-validation procedure were averaged and used for statistics calculations. The performance of models was defined by standard evaluation metrics R 2 and RMSE. Outliers are defined as points that have error in prediction that is more than 3*RMSE of model. This chapter is devoted to the modeling of the logarithm of reaction rate (logk) of bimolecular nucleophilic substitution reaction. Bimolecular nucleophilic substitution reaction is very common and experimentally well-studied type of reaction. Nucleophilic substitution (SN) is a fundamental class of reactions in which an electron rich molecule called nucleophile attacks the positive or partially positive charged atom of substrate molecule to replace a leaving group [197] (see Figure 33.). Bimolecular nucleophilic substitution SN2 is referred to subclass of SN reactions where the bond with leaving group is broken and the bond with nucleophile is formed synchronously. It is important to note that synchronicity is essential feature of SN2 reactions. Nucleophilic substitution reactions that proceed monomolecularly through formation of carbocation with following ion recombination are usually denoted as SN1.

Nucleophiles could be either neutral (usually amine or alcohol) or negatively charged species (alcoholates, thiolates, halogen or other inorganic salt anions, neutral amines). Usually, only reactions with aliphatic carbon in reaction center are called as SN2 reactions (see example on Figure 33.). Reactions that involve substitution at aromatic or unsaturated carbon of substrate atom are usually asynchronous and proceed through addition-elimination (SNAr) or elimination-addition (SN1 or benzyne) mechanism and thus are not called SN2. One should notice that SN1 (two-stage mechanism with monomolecular kinetic equation) and SN2 (one-stage, synchronous with bimolecular kinetic equation) mechanism are indeed two extremes that rarely take place in pure state.

Usually mechanism of reaction is rather complex, having features of both extreme reaction types [198]. Several attempts were made to build models for SN2 reaction rate constant using QSRR approach [START_REF] Varnek | Substructural fragments: an universal language to encode reactions, molecular and supramolecular structures[END_REF][START_REF] Kravtsov | Prediction of Rate Constants of SN2 Reactions by the Multicomponent QSPR Method[END_REF]139,141] In this chapter we report the information about the dataset containing some 8000 SN2 reactions proceeding in 43 different solvents and water-organic mixtures at different temperatures. The data analysis and visualization were performed with the help of Generative Topographic Mapping [199]. For the first time, the Matched Molecular Pairs approach [163] was applied to the analysis of substituent effect. A new technique of for unbiased validation of the structure-reactivity model was suggested. Finally, results of external validation of the proposed model will be discussed.

Data description

Dataset was curated by means of strategy described in section 4.

Data visualization and analysis

Visualization of chemical space is a powerful tool for data analysis of dataset.

However, for chemical reactions this technique was extremely rarely been used since the work of Gasteiger who used Kohonen self-organized maps (SOM) for reaction classification [200]. Visualization by SOM dataset of chemical reactions of different types was used also for selection of reaction signatures [104].

Generative Topographic Mapping (GTM) [201] implemented in our laboratory

[202] have shown its superiority over other visualization methods [199] and have never been used for the analysis of reaction space. In this method, the data points originally located in D-dimensional space (where D is equal to number of descriptors) are projected on 2D latent space (called manifold). The main difference from the other visualization method is that objects are projected probabilistically (with different probabilities, called responsibilities) to grid nodes of manifold, thus every projected object is indeed represented by a distribution on the map. In this case, a position of an object on the map are calculated as gravity centers of its probability distributions. Thus, although manifold is represented by a set of the nodes (grid of points, like in SOM), the positions of object on the GTM map is continuous and not tight to the node positions contrary to SOM.

The power of maps in data analysis relies on the possibility to color objects according to different criteria. Here, we colored data points according to reaction signatures, substrate, and nucleophile nature (see Figure 35). GTM was built using 25x25 grid, with other parameters set to defaults (12 RBF with width 2.8 and regularization set to1) that according to our experience is the best choice. Since we wanted to analyze only structural diversity in reaction space, condition descriptors was omitted, sequences of length from 2 to 4 containing at least one dynamic atom or bond were used as structural descriptors. It should be noted that the reaction rate depends not only on reagents structure but also on experimental conditions. It is clear that MRP reflects only structural factor.

Therefore, only reactions in almost same conditions could be analyzed with MMR. Here, reactions in pure methanol running under ambient conditions (25-35°C) have been selected for this analysis.

For SN2 reaction, substituents effects could be interpreted in the framework of the reaction mechanism where an atom of nucleophile possessing lone electron pair or bearing negative charge attacks partially positively charged carbon atom which results in a leaving group replacement. Thus, electron donating substituents in nucleophile increase its reactivity and, hence, a reaction rate. Similarly act electron-acceptor substituents in substrate molecule which increase partial positive charge on reacting carbon atom. Figure 37 provides several examples of MRP that fully supports the known mechanism and effect of substituents. As expected, substitution of hydrogen in nucleophile by electron acceptor chlorine atom slows down reaction (Figure 37a), replacement of electron donating methoxy-to acceptor nitro-group in substrate molecule increases its speed (Figure 37b). However, some replacements, e.g., hydrogen to methyl group (Figure 37c) lead to small alteration of rate constant in different directions -the sign could be either negative or positive depending on position (ortho-or para-).

Considering that reaction rate is measured with error of almost 0.5 log units (section 4. 1.3.4.), these small logk variations could be attributed to data noise. MRP could also be very helpful for analysis of the data quality. If change of substituent leads to different in sign changes of reaction rate constant that are greater than experimental errors as shown on Figure 37d it could be an indication of plausible error in data. For example, strange trend in logk variation for the MRP describing substitution of hydrogen by nitro-group in substrate (Figure 37d) was observed. One can see that in this case reaction rate constant could either increase or decrease. Indeed, right reaction pair on Figure 37d shows acceleration due to nitro-substituent in substrate molecule and this fact is fully in line with theoretical concepts of substituent effect for SN2 reactions. Moreover, these two reactions have halo-acetophenone substrate for which SN2 reaction mechanism is virtually the most probable. However, in 8 reaction pairs, one of which is shown in left part of Figure 37d, substitution of hydrogen by nitro-group leads to decrease in reaction rate constant. It is common for SN1 reactions which form carbocationic intermediate that is destabilized by electron-withdrawing substituents. However it could also happen to an SN2 reaction with late transition state and great charge separation where the bond with leaving group is strongly loosened. Then partial positive charge on carbon could be destabilized by electron acceptor and thus even in case of SN2 reaction electron withdrawing group could slow down reaction. The primary sources from which these reactions were taken from were carefully examined and we came to conclusion that almost all measurements correspond to SN2 mechanism. However one paper [203] was very old (1925), that time mechanism of nucleophilic substitution reactions were unknown and it seems that unimolecular reaction have been considered bimolecular and thus data were incorrect (strong dependence of rate constant on reagent concentration was found for neutral nucleophile that should not take place if reaction rate constant of SN2 reaction correctly determined). These data points were excluded from the dataset.

Model building and validation

The collected reactions were encoded by CGRs for which ISIDA fragment descriptors were generated. Each solvent was represented by 14 physico-chemical parameters accounted for polarity, polarizability, H-acidity and H-donor ability as well as molar percent of organic solvent in water to model solvent mixtures (it equal to 100% if solvent is pure). Inverse temperature was also used as descrriptor. Variety of fragment descriptors of different size and topology were generated, then concatenated with conditions descriptors and use in the building of Support Vector Regression models.

Optimal SVR hyperparameters and the best fragmentation schemes were selected by genetic algorithm. Ten best fragmentations that allow creation of models with highest predictive performance were selected. For each of them 10 repetitions of 5-fold cross validation have been performed, then all models were saved and used in consensus Analysis of outliers shows that 105 data points are predicted with error more than 3•RMSE. Among those 53 data points correspond to reactions involving anionic nucleophile and 52 to neutral nucleophile.

Further inspection of cross-validation procedure showed that model performance estimation is too optimistic. This situation arises from "naïve" cross validation procedure.

The problem is very straightforward: if two data points correspond to the same reaction proceeding under slightly different conditions, the difference in logk value is small. If one of these reactions is selected to test set the other to training set, the object from test set will be predicted very close to the true value. Hence more similar conditions per transformation are reported, more chance to observe too optimistic estimation of the model performance in cross-validation.

In order to avoid this problem we decided to assess model performance only on data points which were measured under one condition only (we called them unique data points, UDP). As an unbiased estimation of predictive performance one could use prediction of UDP in cross-validation. In this case overestimation of predictive performance is impossible, since a given reaction in the test set can never occur in training set. Among selected 551 UDP, 202 and 349 belong to reactions with anion and neutral nucleophiles, correspondingly. Statistical performance of parameters for UDP reactions in cross-validation procedure is much closer to experimental error: RMSEUDP=0.61 and R 2 UDP=0.75 (see Figure 39 for predicted vs experimental plot on UDPs). Figure 39. Predicted vs experimental values of logk for unique data points. Solid line corresponds to perfect predictions, dashed lines specify margin with values predicted within 3RMSE, crosses and circles are neutral and anionic nucleophile reactions respectively.

Outliers analysis

Examination of data points for which difference between predicted and experimental values exceeds 3•RMSEUDP reveals 18 outliers, among which only 7

reactions with neutral nucleophile and 11 with anionic nucleophile.

Analysis of outliers shows that most of them result from dataset imperfectness and modeling procedure used. The errors were caused by the following reasons (reactions are drawn in Table 6).

• Non-continuous dependency of the rate constant on temperature, reaction 1. E.g. logk for a given temperature -4.09 (90 °C), -3.81 (100 °C), -3.68 (105 °C), -3.58

(110 °C), -3.42 (115 °C). Such a small rate constant -3.21 (201 °C) is out of trend and model logically predict it as -1.36.

• Complex structural effects that were not learned by model due to lack of representatives. Reaction 5 represents effect of direct polar conjugation for nucleophile, which strongly influences reaction rate however could not be captured by a model since it is the only reaction with para-nitroaniline.

• Transformations possessing rarely occurred structural patterns in dataset are outliers (items 16, 11, 6, 3), since the number of instances possessing given fragment in the training set is not sufficient to learn the contribution of these patterns. E.g. there was only 1 substrate possessing -SO2CH2Cl group (item 16)

thus it was mispredicted. Such a large rate constant for reaction 11 could be explained by anchimeric assistance, due to lack of data this effect is fully neglected by the model. Reaction 3 has rare transformation of tertiary amine to quaternary.

The opposite effect could happen for small size reactants constituted of common small fragments (reactions 2, 8, 12, 17, 18). These fragments can hardly distinguish the difference between the reactant structure.

• Solvent could have unexpected influence on the rate constant that is not captured by the model. For example, transformation 7 has 3 magnitudes higher logk in aprotonic polar DMFA solvent than in protonic methanol and ethanol solvents present on training set. Non-continuous dependency of the rate constant on solvent mixture composition leads to outlying prediction for reaction 10.

According to training set data rate speed increases with rise of DMSO percentage in mixture with water. For instance,. logk for a given percentage of DMSO in water (given in brackets) is the following: 1.07 (81% DMSO), 0.06 (65% DMSO), -1.01 (46% DMSO), -2.08 (30% DMSO), -2.93 (18% DMSO), -3.66 (8% DMSO), -4.09 (2% DMSO). Thus, such a small rate constant for pure DMSO -2.48 (100% DMSO) is out of the trend.

• Mistakes in data source. In primary source (paper [204]) rate constants for items 4, 9, 13 correspond to reactions of benzylamines while in reference book [START_REF] Palm | Tables of Rate and Equilibrium Constants of Heterolytic Organic Reactions[END_REF] from which data were taken reactions were annotated as phenylamines. Thus in our database wrong substrates were annotated. Alkylamine group is stronger nucleophile than amine conjugated with phenyl ring and the reaction rate constants for the reaction with the latter about 3 orders of magnitude smaller.

• Another problem revealed from the outlier analysis is based on reactions with similar descriptors but different properties. If two very similar reactions with drastically different logk are present in the dataset, the both can be mispredicted. The question arises: whether solvent descriptors correctly capture solvent effects ?

+ + + + H 3 C O Cl O O O H 3 C N CH 3 H 3 C CH 3 H 3 C N + H 3 C CH 3 O - Cl O O O + + CH 3 I N H 3 C CH 3 Cl CH 3 N + CH 3 H 3 C Cl I - + + H 3 C O Cl O O O -O NH 2 O N + O - H 3 C NH O N + H O Cl O O O + + Br O O O H 3 C H 3 C H 3 C N O O O CH 3 CH 3 CH 3 N + Br - + + CH 3 I O - Na + H 3 C O I - Na + + + + + H 3 C Br N + O O - -O Ag + CH 3 N + O O -O
From the entire set, we selected 6 subsets of reasonable size for the reactions proceeding in particular solvents (nitrobenzene, methanol, ethanol, acetone, water, benzene). The models were built on these subsets using only fragment descriptors and unique data points. In most of cases, RMSE obtained in cross-validation for particular datasets is similar for the global and local models (Figure 40). This means that solvent descriptors are rather good to account for solvent effect in logk modeling . Notice that accuracy of predictions is not similar for different solvents: prediction error observed for nitrobenzene and ethanol subsets is smaller than for other solvents

Validation on the external set

For validation of the model, external data set containing 104 Menshutkin reactions was collected from the papers published in 1990-2010s. Since the reference book [START_REF] Palm | Tables of Rate and Equilibrium Constants of Heterolytic Organic Reactions[END_REF] serving the data source was dated 1978, external set data don't overlap with the training set ones. Predicton performance was slightly worse than that observed in cross validation for UDP (see Figure 41) RMSE=0.8 and R 2 =0.64. We tried several applicability domain definitions but none of them was perfect. Thus, selecting reactions for which reaction center with its second environment equivalent to that in training set reactions (similar to ICClassify narrow signatures [205]), Figure 41b, retains only 13% reactions with high RMSE = 0.69 and low R 2 = 0.22. At the same time reaction center signature with first environment considers all reactions lying within AD.

The "consensus control" applicability domain [160] was also considered.

According to its definition, consensus prediction is considered unreliable if a given reaction is outside of bounding box AD for a certain percentage (50% by default) of individual models. This AD was too restrictive although it efficiently discards the outliers.

If we smooth the consensus control requirements and accept even one individual model considering a given reaction within its AD, this leads to RMSE = 0.61 and R 2 = 0.5 with 32.3% coverage (Figure 41c). Raising threshold to 50% leads to retaining very few similar to training set reactions for which, however, logk was perfectly predicted (RMSE = 0.1, R 2 = 0.98, coverage = 7.5%), Figure 41d.

Conclusions

The consensus model for the rate constant of SN2 reaction proceeding under different reaction conditions has been built using fragment descriptors generated for Condensed Graph of Reactions and special descriptors accounting for experimental conditions. The model displays a reasonable performance both in cross-validation and on the external test set. The global model obtained on the entire set performs similarly to local models built on the subsets corresponding to particular solvents or nucleophile types. The models are available for the users on our server (cimm.kpfu.ru).

We have demonstrated that Matched Reaction Pairs approach could efficiently be applied for the analysis of substituent effect. It was found that mostly it is fully in line with theoretical concepts issued from the reaction mechanism. Detection of unusual substituent effect is in MRP analysis, could help to identify either unusual reaction mechanism or error in data annotation. Thus, both data visualization and MRP analysis could be used as tools facilitating data cleaning process. E2 reactions were used in the benchmarking study. This dataset is not as large as SN2 one and reliable results could be obtained in reasonable time. Some of descriptor spaces (for example, for based on SiRMS descriptors of mixtures) were extremely large and descriptor storage was also an important issue. Moreover, it was decided to reduce number of reactions in data set by careful manual examination of the E2 dataset and exclude all doubtful data, data where stereoisomery was important, data containing structural errors. This was done to ensure that the descriptors indeed reflect the relevant The modeling set contained 880 Diels-Alder reactions measured only in pure solvents. In cycloaddition not only initial stereo configuration of reagent molecules, but also endo and exo cycloaddition are characterized by different rate constants. So, even one reaction with specified stereo configuration of reagents can give two products.

Typically this information was neglected in the initial sources, but for some reaction stereoconfiguration of reagents and products was known. For the latter we found the difference between rate constants values, which, fortunately, didn't exceed interlaboratory error of some 0.5 logk units. That's why, the data on stereoisomers were merged and rate constants were averaged. The modeling procedure and results are described in the article published in Russian Journal of Structural Chemistry, see below.

Chapter 8.

Modeling of tautomeric equilibrium constants

Tautomerism represents one of the most common phenomenon in organic chemistry. According to IUPAC definition, the tautomerism is the easily undergoing isomerism of the following general type

The most common form of tautomerism is prototropy (prototropic tautomerism),

where G is hydrogen. The importance of tautomerism in drug design was pointed out in numerous publications [129,[207][208][209][210][211][212]. The main problem with tautomers is related to the fact that each tautomer behaves as an individual chemical compound with its own spectral, chemical and physical properties and, at the same time, all tautomers of a given compound coexist in dynamic equilibrium. For this reason, identification of molecules and its properties that are able to tautomerize can be challenging.

Since a molecule can exist in several forms with different descriptor vectors and structure fingerprints, one have to decide which form should be taken into account in the course of modeling or storage in chemical database. It has been reported several times that several tautomeric forms of the same compound can be found in the same database, and even different prices are given for them [213,214]. In QSAR studies and chemical database management only the most stable (most populated) tautomeric form in water is considered mainly. It should however be taken into account that guessing the most stable form may sometimes be rather difficult problem, because tautomer equilibrium is known to be solvent and temperature dependent. Selecting one of several alternative tautomers is usually performed using rules of tautomer canonicalization [215][216][217]. In structure-based drug design several tautomeric forms with energy lying within the given energetic window are considered. In this case relative stabilities of tautomers can be estimated using quantum-chemical [218,219] or force-field modeling [220,221].

The dataset of tautomeric equilibria was collected from reference book [START_REF] Palm | Tables of Rate and Equilibrium Constants of Heterolytic Organic Reactions[END_REF]. The tautomeric equilibria were represented as reactions that turn one compound to another, in order to fit equilibria to the developed approach. In this representation equilibria constants belong to one exact direction of process. This required to check additionally that left and right parts of the reaction equation are represented in a unified way. E.g.

keto-enol tautomerism have to be represented only as transformation of ketone into enol.

Additional problem we came across is, that using dynamic bond concept in initial

X Y Z G Y Z X G
formulation of CGR [START_REF] Varnek | Substructural fragments: an universal language to encode reactions, molecular and supramolecular structures[END_REF] with implicit hydrogens, encoding of neutral -zwitter-ion tautomeric equilibria produces CGR without any dynamic bond. Therefore, dynamic atoms were introduced for description of such reactions.

Initial curated set of 840 records was curated for modeling, duplicated reactions' rates were averaged and resulted modeling set was divided into training and external set. Taran contained 18 reaction. In this chapter we will describe the problems we faced and the solution proposed. [6,235] as well as alkyne nature [233] were shown to have great impact on the reaction rate switching it from moderate to ultra-fast. Very recently detailed stopped-flow kinetics study of reaction of some sydnones have been performed that supported twostage mechanism and rates of cycloaddition and retro-Diels-Alder steps were measured [6]. It was shown that for fast reactions of fluorosydnones the second, CO2 release step, was rate-limiting. Thus, by now, major attention was devoted to design of alkynes or alkenes with improved reactivity toward strain promoted 1,3-cycloaddition. In this part of the work we wanted to create predictive model for reaction rate of sydnones with alkynes and to shed light how substituents in sydnone ring influence its reactivity. We collected some experimantal data on sydnone-alkyne cycloaddition that was used for assessment of the error of computational prediction of reaction rate and revealing the importance of electronic and steric effects of substituents of sydnone ring for reaction rate.

Data set description

The data on cycloaddition reaction rate were taken from publications of F. Taran's group [229,234], Figure 44. The collected dataset contains reaction rates at ambient temperatures of 15 reactions of sydnones and BCN, 3 reactions of sydnones with TMTH (see Table 7).

QSRR modeling

In the beginning of the project the attempt to model rate constants using QSAR approach was done. To do it some new, confidential data on reaction rate with BCN was given to us by F. Taran's group. Totally we had 34 data on measurement of rate constants of different sydnones and iminosydnones with BCN. The only factor influencing the reaction rate was the chemical structure of the sydnones, since the alkyne and reaction conditions were kept constant. Thus, there is no need in application of CGR approach for the reaction rate prediction and standard QSPR approach could be used. However all approaches have shown poor performance (Q 2 <0.5) on 10-fold crossvalidation. We explain it by following factors:

• Small and heterogeneous data set -the number of points in the dataset is not enough to learn all the factors that influence on reaction. The data set appeared to be too heterogeneous, most reactions are slow and 4 had a much bigger reaction rate constant than most of the others in the dataset.

• Obvious feature that influence reaction speed -one can notice that reaction rate increases when halogen atom is present in X position. Thus upon descriptor selection usually presence of certain halogen becomes an important feature.

However such model could not predict correctly test set reactions where a sydnone contained another halogen atom. Such effect and small number of high-speed representatives makes impossible to learn other features from this dataset.

Our attempt to manually select descriptors according to knowledge of mechanism failed as well and we never managed to build a model with moderate performance. Thus we came to conclusion that the dataset has low modelability for QSPR approach application due to its composition. Thus, we turned toward quantum-chemical calculations to estimate the reaction rate.

The work was divided into two parts: development of the workflow for direct calculation of transition state Gibbs free energy with following reaction rate constant calculation and revealing the structural effects responsible for reaction rate. The tasks were solved using dataset given in Table 7.

Details of quantum chemical calculations

In this sections the computational approaches used for calculation will be described.

Energy and geometry optimization

Density Functional Theory (DFT) [118,119]

𝜔 = 𝜇 = 2𝜂
where μ is electronic chemical potential, and η -chemical hardness. They could be expressed in terms of HOMO (εHOMO) and LUMO (εLUMO) energies as:

𝜇 = (𝜀 ©ª«ª + 𝜀 ¬-«ª )/2 𝜂 = 𝜀 ¬-«ª -𝜀 ©ª«ª
The HOMO and LUMO energies were obtained within DFT scheme [118] for sydnone molecules.

Reaction pathway investigation

The full reaction path was explored for reactions 2, 15 and 16 in Table 7 whose kinetic measurements using stopped-flow technique was performed in reference [6]. The One can see from Figure 47 that the limiting step of reaction 15 is the first transition state (TS1), intermediate is lower in free energy than reagents by 21 kcal/mol. According to our calculation, intermediate is very unstable since it is separated from products by transition state (TS2) with a tiny barrier (some 1 kcal/mol). Intermediate decomposition is exergonic by some 85 kcal/mol. The same is true for the much slower reaction 2.

Solvent effect that was accounted by IEF-PCM model [START_REF] Tomasi | The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level[END_REF] shown almost no influence on the reaction barrier due to compensation. In reaction 16 between fluorosydnone and TMTH TS2 was not localized at all. Liu at al [6] experimentally observed the existence of an instable intermediate during the reaction. However, they estimated the second step reaction rate in the range 0.02 to 0.98 s -1 depending on alkyne and concluded that this second step becomes rate limiting. The discrepancy with the computed free energy clearly points toward yet undefined entropic effects that are balancing the two steps of the reaction. A possible hypothesis is that the solvent water molecules play an active role in the studied reactions.

With the present set of approximations, our results show that there is no qualitative difference between reaction paths for fast reaction of fluorosydnones (X=F) and much slower reaction of unsubstituted (X=H) sydnones: both are characterized by an unstable intermediate. This is also in qualitative agreements with the observations reported by Narayanam et al [233] for unsubstituted sydnones with strained alkynes.

Assessment of rate constant

The study of Liu et al [6] reported that the first step of the reaction is strongly affected by the nature of the sydnone. The main result is that the rate of the first step of the reaction is so fast, if it involves a fluorosydnones, that the second step became the limiting. This second step seems marginally affected by the nature of the sydnone and much more by the nature of the strained alkyne, ranging from 0.02 s -1 to 1 s -1 . However, these variations are order of magnitudes smaller than those that affect the reaction rate of the first step when varying the nature of the sydnone, ranging from 42 Mole -1 .s -1 to 12000

Mole -1 .s -1 .

This explains the paradoxe that the effective rate of reaction (rate of product formation) seems to be mostly driven by the nature of the sydnone and the rate of the first step (see Table 7). The first step is so fast and the transformation is so irreversible that it produces a large accumulation of the instable intermediate. In turn, this produces a massive imbalance that pushes toward the generation of the final products. Moreover, according to description of HPLC experiment in author's previous studies [235] rate constants reported in Table 7 represent the speed of conversion of reagents. The latter depends only on first step rate constant.

Hence, from practical point of view the first step rate constant is the most relevant.

Experimental measurements of rate constants are time-consuming and expensive while the variability of structure is rather big. For practical application it is important to find a way to quantitatively assess the rate constants of some reaction. We attempted to find activation free energies for all reactions under study and compare predicted values with experimental ones. The goal was to assess the quality of quantum-chemical estimation of reaction kinetic characteristics. For all reaction under study TS1 was identified and its activation energy was calculated as difference between free energy of TS1 and reagents. However, due to large number of reactions under study, conformational lability of interacting compounds and many possible orientations of reagents to form transition state, an automatic approach for the detection of TS1 was needed (described in details in Chapter 4. Free energies of activation, ΔΔG ≠ calc, were predicted for all reactions under study using described quantum chemical approach. For comparison, "experimental" values of free energies, ΔΔG ≠ exp, of activation were calculated using Transition State theory and Eyring equation (see Chapter 4.3.3.3) from measured rate constants. The comparison of these two quantities is given on Figure 48. We suggested that sydnone-alkyne reaction could be charge controlled. Fukui nucleophilicity, electrophilicity, radical attack susceptibility using Conceptual DFT [244] for 5 core atoms of sydnones were calculated and used to find correlations with logk of sydnones -BCN cycloaddition. Only one unbiased correlation was found between logk and nucleophilicity of C3 atom with R 2 = 0.43 (Figure 50b). However, one could notice that correlation was mainly caused by fluorosydnone-BCN reaction 15. The deletion of the latter point will lead to dramatic drop of correlation coefficient. Thus, charge control does not play the major role in these cycloaddition reactions.

The result leads us to conclusion that there is a complex interplay of structural factors that cannot be caught by simple linear correlations. To avoid any additional effect reaction rate in small congeneric series of reactions of sydnones having different substituents X and the same R with BCN (reactions 3, 10-12, 15) and fluorosydnones (X=F) having different R with TMTH (reactions 16-18) was analyzed.

Results for reaction of BCN with halo-sydnones [START_REF] Drucker | Support vector regression machines[END_REF][START_REF] Blakemore | Organic synthesis provides opportunities to transform drug discovery[END_REF][START_REF] Peverati | Quest for a universal density functional: the accuracy of density functionals across a broad spectrum of databases in chemistry and physics[END_REF][START_REF] Cohen | Quantities, Units and Symbols in Physical Chemistry[END_REF][START_REF] Manion | NIST Chemical Kinetics Database[END_REF] given on Figure 51 in blue show significant variation of electrostatic charges on С3 atom, which change from negative (X = H) to positive value (X = F), whereas the charge on N2 atom slightly varies as a function of substituents. In this series, the logarithm of the cycloaddition rate (logk) well correlates with the C3 charges (Figure 51a), r=0.98 and RMSE = 0.13 logk, p value 4•10 -4 . At the same time no correlation with LUMO energy was observed (Figure 51b).

Variation of aryl substituents at N3 atom in sydnones in reaction 16-18 (with TMTH, in red on Figure 51a), doesn't lead to significant variations of charge at C3 atom (Figure 51a). One could notice that the rate of reaction with TMTH is affected to charge variation in greater extent than BCN. Thus electron nature of substituents on sydnone ring will have stronger influence on reaction with TMTH. On the other hand, LUMO energy in this series decreases from -3.57 to -3.05 eV which could partially explain considerable raise of the rate constant (Figure 51b). One should take into account that correlation coefficient is not big enough (0.86) and 3 points are too few to make robust conclusions. 6, 353.2, 351.9, 351.5, 351.2 correspondingly. This measure can evaluate steric tension in TS1. As one can see sydnone with F is even more planar in TS1 than sydnone with H. We can consider that negative charge on F atom and positive charge on hydrogens of alkyne can form weak hydrogen bond in this phase, that reduces energy of TS1.

These observations are opposite to suggestions by Liang et al [236] that electronically more reactive electrophiles should be sterically more encumbered. Indeed, on one hand, iodo-sydnone, sterically more encumbered than fluoro-sydnone, is less reactive than the latter due to weaker electron acceptor ability and, on the other hand, steric repulsion of halogen with alkyne molecule. But, indeed, fluorine change to more bulky electron-withdrawing substituent could lead to drastic loss of reaction rate. Thus, reaction 9 involving sydnone with X=CF3 is drastically slower than 10 with X=Cl, despite atomic charges on C3 are close (C3 charge is 0.0093 when X=CF3, 0.0042 for X=Cl, 0.0891 for X=F).

In the context of rational design of highly reactive substrates, our calculations result in the following conclusions:

1.

Strong electron acceptor X at C3 atom increases its positive charge and, thus, improves its affinity to alkyne in cycloaddition. However, this substituent should not be bulky because it may cause steric hindrance to alkyne in the transition state. In this context, fluorine substituent at C4 represents an optimal choice. However, other flat acceptors could be a good alternative as well.

2.

Strong electron acceptor at N3 atom weakly affects charge distribution in sydnone moiety, but decreases its LUMO energy, which, in turn, favors cycloaddition, especially in reaction with TMTH. Thus, one can expect that nitro-groups or other strong electron withdrawing substituents in benzene ring R will strongly favor reaction.

3.

Compared to BCN, more sterically strained TMTH has higher HOMO energy which explains its better reactivity. the same systems by Liu et al [6]. They discovered that the second step of the reaction was the limiting step, in contrast to our results. Yet, for practical application, the effective reaction rates are reporting the consumption of the reagents which is controlled by the rate of the first step of the reaction. Besides, this rate of the first step is much more variable with the nature of the reagents than the second step. This is confirmed by our quantum-chemical calculations well reproducing experimental data on reaction kinetics based on Transition State theory. The root-mean squared deviation between quantumchemically calculated activation free energies and the value that was recalculated from experimental data on reaction rate constant using Eyring equantion on 18 reaction of sydnone-alkyne cycloaddition is lower than 2 kcal/mol. We did not have yet a satisfactory understanding of the mechanism of the second step of the reaction. But, since our results are not changed using an implicit solvation model, it is possible that some solvent water molecules play an active role in the process.

From the other hand, such a good reproducibility of activation free energy and as a consequence reasonable prediction of rate constant by quantum chemical calculations gives a powerful tool for relatively cheap screening of possible reagents. Application of developed workflow could be used to computationally prove hypothesis on reactivity of certain pairs of sydnone and alkyne. Only the most promising candidates could be synthesized and their reaction rate could be experimentally measured.

To reduce search space in this work the structural factors important for reaction rate were analyzed. The absence of clear correlation within whole dataset of reactions could be explained by complex interplay of effects responsible for reaction rate and some experimental noise in reaction rate constant measurement. However, in restricted series of reactions some correlations could be found that made possible to reveal three major factors affecting reaction rate. First, large positive charge on C3 atom guarantees faster reactions (atom numbering according to Figure 47). This charge itself is mainly affected by electron withdrawing ability of substituent at this atom (X) and partially on the electron withdrawing ability of substituent at N2 atom (R). On the other hand, bulky substituents at C3 atom lead to sterical hindrance to approaching reagents and complicate transition state formation. Thus, only sterically unencumbered electron withdrawing groups could be used at C3 atom to gain reaction speed. Since LUMO of sydnone is participating in orbital interactions, its energy is an important factor influencing rate constant. The lower the energy of LUMO the greater is rate of reaction.

Mostly LUMO energy is affected by π-electron withdrawing ability of substituent R.

Thus, usage of stronger electron acceptors will favor reaction. Being augmented by proposed workflow for fast screening of reaction rate these recommendations provide an efficient tool to the design of more active agents for bioorthogonal click reactions.

The quality of correlation between free energy of activation calculated quantumchemically and estimated from experimental value of rate constant using Eyring equation

shows from one side that developed quantum-chemical approach reproduces the energy rather well in absolute scale (less than 2 kcal/mol). From the other side, having such good description the correlation coefficient is not large (r is 0.84, that means that determination coefficient would be at most 0.6). The reason could be noise in data and in this case it will explain why we failed to build QSPR model on the dataset -high level of noise in addition to other problem such as imbalances and heterogeneity of dataset prevented building predictive model.

1. The first database of chemical reaction kinetic and thermodynamic properties was collected. Information on more than 10 000 chemical reactions involving structure, solvent, temperature, and rate constants of bimolecular nucleophilic substitution, bimolecular elimination, cycloaddition (Diels-Alder reactions) and tautomeric equilibrium constants was annotated from reference books and PhD thesis defended in Kazan Federal University.

2. The workflow for reaction standardization and curation was proposed and the tools required for it was developed. Some tools are based on approach of Condensed Graph of Reaction. The latter was extended by notion of dynamic atom to extend its applicability to wider range of reactions. The approach to store CGRs in SDF format was proposed. The toolbox developed includes CGR hashing, reaction center detection and hashing, workflows for AAM checking and conditions verification. 6. Predictive models for rate constants of bimolecular nucleophilic substitution, bimolecular elimination, and Diels-Alder reactions, as well as tautomeric equilibrium constants were built. The models were published on-line on server developed specially for chemical reactions.

7. To predict sydnone-alkyne cycloaddition reaction rate constant the workflow based on quantum-chemical calculations and semi-automatic identification of transition state was developed. QSPR modeling of this dataset failed. Using quantum chemistry approach, activation free energy of reactions under study were reproduced with some 2 kcal/mol accuracy. To speed up selection of optimal reagents for these cycloadditions, most important factors affecting the reaction rates were reported.

  SIlico Design and data Analysis SMIRKS-Linear representation for reactions RInChI -International Chemical Identifier for Reactions (alphanumerical string for reaction structure representation) LFER -Linear Free Energy Relation IUPAC -International Union of Pure and Applied Chemistry SMILES -simplified molecular-input line-entry system DFT -Density Functional Theory HPLC -High-performance liquid chromatography SiRMS -SImplex Representation of Molecular Structure MRP -Matched Reaction Pair KFU -Kazan Federal University LS -Latent Space CRDB -Comprehensive Réaction Database SN2 -Bimolecular Nucleophilic Substitution reactions E2 -Bimolecular Elimination reactions DA -Diels-Alder reactions

Figure 1 .

 1 Figure 1. Exemple d'une réaction chimique (à gauche) et du Graphe Condensé associé (à droite). Les labels et correspondent à des liaisons créées et cassées, respectivement. Les labels "C-1" et "C+1" décrivent les charges dynamiques caractérisant les atomes dont la catégorie de charge (neutre, positive, négative) varie au cours de la transformation chimique.

Figure 2 .

 2 Figure 2. Carte GTM sur le jeu de données de 1394 transformations de réactions SN2, encodées par des fragments ISIDA, sans prendre en compte les descripteurs de conditions. Les objets sont colorés en fonction a) de la signature du

Figure 4 .

 4 Figure 4. Corrélation entre les énergies libres d'activation calculées (Predicted ∆∆G) et les observations expérimentales (Experimental ∆∆Gcalc) pour les réactions entre les sydnones et les iminosydnones avec des alcynes cycliques. Les ∆∆G expérimentales sont calculées à partir des valeurs mesurées de logk en utilisant l'équation d'Arrhenius. Les coefficients de détermination (R 2 ), les erreurs quadratiques moyennes (RMSE) et le nombre de points de données (n) sont affiché dans le graphique.

  Figure 5) allows the encoding of a reaction as pseudo molecule and to apply all QSPR techniques for modeling of its properties.The thesis consists of 10 chapters. Chapter 2 provides information about the approaches and tools of chemoinformatics and quantum chemistry used in this study.The 3 rd chapter describes different types of reactions encoding in chemoinformatics and reviews publications on structure-reactivity modeling and provides with some information about chemoinformatics and quantum chemistry approaches and tools used in this study.[START_REF] Tomasi | The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level[END_REF] th chapter describes the Comprehensive Reaction Database (CRDB) created in this study and procedures of data curation, cleaning and new way of representation of reactions with the help of CGR signatures and dynamic atoms. The chapter 5 describes predictive models for the rate constant for SN2 reactions. The Next chapter 6 is devoted to modeling of E2 reactions and benchmarking of ISIDA descriptors in comparison with classical ones and new type of mixture SIRMS descriptors for reactions. In the 7 th chapter application of the approach to cycloaddition reactions in solution is described. Chapter 8 is devoted to the SVR and quantum chemical modeling (QC) of the tautomeric equilibrium constants. In Chapter 9, the SVR and QC modeling of the rate constant of several bioorthogonal reactions involving sydnones is described. The Chapter 10 summarize conclusions for the work.

Figure 5 .

 5 Figure 5. Example of a chemical reaction (left) and related Condensed Graph (right). Labels "0>1" and "1>0" correspond to created and broken bonds, respectively.

[C: 1 ]

 1 (=[O:2])[Cl:3].[H:99][N:4]([H:100])[C:0]>>[C:1](=[O:2])[N:4]([H:100])[C:0].[Cl:3][H:99]

Figure 7 .

 7 Figure 7. Reaction (top) and corresponding SMIRKS (bottom). Atom-to-atom mapping is specified by numbers.

Figure 8 .

 8 Figure 8. RInChI for a given concurrent reaction.

Figure 9 .

 9 Figure 9. Common representation of chemical reaction with SN2 mechanism.

Figure 10 .

 10 Figure 10. Comparison of different reaction representations.

Figure 11 .

 11 Figure 11. Common representation and CGR representation of reaction with addition features (dynamic bonds and dynamic charges). Bond with circle -dynamic bond corresponding to formed single bond, crossed bond represent cleaved single bond. "c+1" and "c-1" labels are used to represent dynamic atoms: increase or decrease of atomic formal charge by one, respectively

Figure 12 .

 12 Figure 12. Ugi-Dugundji R matrix as representation of reaction.

Figure 13 .

 13 Figure 13. Vector of numbers after ISIDA fragmentor.

Figure 14 .

 14 Figure 14. Generation of SiRMS descriptors for chemical compound.

Figure 15 .

 15 Figure 15. Example of fragment descriptors generation for CGR. One could notice that the procedure is very similar to the one used for molecules descriptor generation, given on Figure 13.

  polarizability. The dielectric permittivity influence on a dissolved molecule is a highly nonlinear. Reaction field theory (a.k.a Onsager theory) [106] gives some formulae that describe the dependence of solute energy on dielectric permittivity: Born and Kirkwood , Debye function [107]. Analogously polarizability of solvent could be described as a function of the refractive index measured by D-line of sodium spectra at 20 degrees Celsius [107]. Several physicochemical scales for description of solvent effects (polaritypolarizability, H-acidity and basicity) were proposed: Catalan SPP, SA and SB scales [108-110], Camlet-Taft π * , α and β scales [111-113], Koppel-Palm tetraparametric equations [114] for solvent effect modeling. Most of these scales are based of solvatochromic effects of the solvent, i.e. shifts in UV, visible or IR radiation absorption

  regression and physicochemical descriptors are used almost exclusively. In SAR/QSAR/QSPR black box machine learning methods based on very flexible fitting functions 𝑓(𝑿) and many other type of descriptors are widely used. The main accent is done on general applicability and robustness of the model rather than its interpretability. At the same, the other methods are widely used for prediction of chemical objects 20 D n characteristics, first of all one need to mention quantum chemistry and molecular mechanics. These approaches use strict physical theories to predict desired properties of chemical systems.

  in house software developed by Dr. Pavel Polishchuk available at https://github.com/DrrDom/sirms. The ISIDA substructural fragments based on CGR were commonly used in the work. ISIDA Fragmentor tool with the following settings were used for the descriptor set generation (see chapter 2.3.1.1 for options description):

  All fragmentations were supplemented by 13 descriptors of solvents. Each solvent was described by special descriptors that represent polarity, polarizability, H-acidity and basicity: Catalan SPP, SA, and SB constants [108-110], Kamlet-Taft constants α, β, and π* [111-113], four functions describing solvent polarity depending on the dielectric constant ε (

  validation procedure. Cross-validation is a model validation technique for assessing predictive power of the model. In 5-fold cross-validation procedure the whole dataset is divided into five almost equal parts. Four of them were used for model building (training set) and one part is used to test predictive performance (test set), see Figure16. The procedure is repeated in such a way that every part once is used as test set. The is done by taking into account only predictions for the reactions when they were in the test set.

Figure 16 .

 16 Figure 16. Cross-validation procedure of machine learning parameters evaluation.

Figure 17 .

 17 Any hyperplane can be written as the set of points 𝑥 satisfying condition 𝒘 • 𝒙 -𝑏 = 0 where 𝒘 is the (not necessarily normalized) normal vector to the hyperplane. The parameter t ||𝒘|| determines the offset of the hyperplane from the origin along the normal vector [152], Figure 17. Maximization of margin was shown to be equivalent to maximization of = ‖𝒘‖ . One could notice that correct classification of objects by plane with given w and b means that 𝑦 j (𝒘 • 𝒙 𝒊 -𝑏) ≥ 1

Figure 17 .

 17 Figure 17. Maximum-margin hyperplane and margins for an SVM trained with samples from two classes. Samples on the margin are called the support vectors.

Figure 18 .

 18 Figure 18. Illustration of parameters of SVR method.

Figure 19 .

 19 Figure 19. Schematic representation of the GTM algorithm: nonlinear mapping of grid nodes in 2D latent space onto a manifold in D-dimensional data space.

  One of the most simple applicability domains is bounding box [158]. It states that an instance is out of model's applicability domain if descriptor values for it are not within the min-max ranges valid for training set. The Bounding Box techniques by definition encompasses so-called Fragment Control: if the data set encoded in structural fragment descriptors, then any molecule of the test set possessing a new structural fragment considered to be out of AD.

Figure 20 .

 20 Figure 20. Representation of the Applicability Domain boundary in chemical space.

Figure 21 .

 21 Figure 21. Example of MMP due to substitution of hydrogen by chlorine.

Figure 22 .

 22 Figure 22. Example of MMP for reactions encoded by CGR.

Figure 23 .

 23 Figure 23. Interactive creation of reaction with Marvin Web application.

Figure 26 .

 26 Figure 26. Dynamic atom representation in tautomeric equilibria. (ChemAxon Marvin Sketch representation [180]). Left side original tautomeric transformation, right side CGR representation with dynamic atoms.

Figure 27 )Figure 27 .Figure 27 .

 272727 Figure 27). Long CGRS text strings could be hashed for shortness and represented as hashed hexadecimal number called CGRS-key (Figure24).

Figure 28 .

 28 Figure 28. Unclear result of reaction in literature.

  ) structure and transformation cleaning, (2) experimental facts cleaning. First step produces standardized transformation representation with mechanistically correct atom-to-atom mapping. The second step is required for discarding duplicates, producing correct tuple transformation-conditionproperty for modeling. Both steps are based on some conventions. For structure cleaning one should decide how the product is specified, since the chemically correct representation is not implicitly best for the modeling purpose. So, it is not really important that structure should be chemically correct. It is more important, that similar transformations have similar representations. Since different people extracted reactions manually, one have to take into account possibility of different representations of molecules participating in reaction, existence AAM errors. These issues should be fixed.For the second step one has to decide what condition parameters influence rate constant. Our convention was that concentration of reagents and additives should not influence the reaction rate constant. And thus reactions proceeding under the same temperature and solvent but with different initial concentration of reagents should be considered duplicates.

Figure 29 .

 29 Figure 29. Conventional reaction representation (left) and CGR (right). Cleaved bonds are crossed, formed bonds are denoted by circle, dynamic atom option c+1 and c-1 mean that charge of atom increased and decreased by one respectively. Numbers near atoms represent AAM. Examples of reaction with formally different but correct AAM (A and B) and reaction with wrong AAM (C) are shown. CGRS and reaction center signatures (RCS) are shown below reactions. For human readability CGRS are compressed using MD5 hash function [194] and represented as hexadecimal number CGRS-key.

Figure 30 .

 30 Figure 30. Difference in logk values for duplicated reactions.

Figure 31 .

 31 Figure 31. Structure of the reaction database. Calculated fields are shown in bold, required fields by plain text and optional fields in italics.

Figure 32 .

 32 Figure 32. General modeling workflow

Figure 33 .

 33 Figure 33. Example of SN2 reaction.

Figure 34 .

 34 Figure 34. Data distribution with respect to (A) rate constant, (B) temperature, (C) solvent (dark part is for pure solvent one and light part is for mixture with water), (D) experimental conditions per transformation.

Figure 35 .

 35 Figure 35. GTM map on 1394 transformations encoded by ISIDA fragments. Objects are colored according to a) reaction center signature (only reaction center atoms included), b) substrates, c) nucleophile structure, d) nucleophile type. The most popular signatures or molecules are explicitly shown.

  known and widely used approach in medicinal chemistry [163]. MMP is defined for a pair of molecules, which are different with a respect of a single group. The extension of MMPs to chemical reactions encoded by CGRs is straightforward since CGR represents a molecular graph with additional atom and bond labels. Thus, instead of comparing a pair of compounds, one can compare a pair of reactions which we'll further call Matched Reaction Pairs (MRP), see Figure 36. It allows to understand how a variation in structure of reactants influences a speed or other property of chemical reaction.

Figure 36 .

 36 Figure 36. MMP built with the OCHEM software for molecules (left) and reactions encoded by CGR (right).

Figure 37 .

 37 Figure 37. Examples of Molecular Reaction Pairs for SN2 reactions in methanol at ambient temperature. In CGRs, formed and broken bond are shown as crossed and circled respectively. a) Reaction rate decreases due to substitution of hydrogen by chloride groups in nucleophile. b) Reaction rate constant increases due to replacement of metoxy to nitro group in substrate. c) Replacement of hydrogen to methyl group leads to small changes in logk depending on the position of the group in aromatic ring of aminoaromatic nucleophile. d) Substitution by nitro-group in substrate leading to decrease (left) or increase (right) of the rate constant.

  Figure 38.

Figure 38 .

 38 Figure 38. Predicted (with global model) vs experimental logk values. Solid line correspond to perfect predictions, doted lines specify margin with values predicted within

Figure 40 .

 40 Figure 40. Cross-validated RMSE of global and local models on the subsets corresponding to particular solvents.

Figure 41 .

 41 Figure 41. External set prediction by the global model. Solid line corresponds to perfect prediction. Dotted line are 3*RMSEUDP away from perfect prediction. (a) reactions are labelled according to signature types. (b) Only reactions having reaction signatures similar to training set ones are shown. (c) only reactions within bounding box AD for, at least, one individual model are retained. (d) only reactions for which 50% of individual models were retained by the bounding box AD are shown

6. 1 .

 1 Models built on CGR-based reaction descriptorsBimolecular elimination reaction (E2) is base-assisted simultaneous reaction of cleavage of bonds with electron-withdrawing group and hydrogen near single bond resulting in formation of a double bond. As for SN2 reaction, it involves a one-step mechanism in which carbon-hydrogen and carbon-leaving group bonds break simultaneously and kinetic equation has second-order (first order with respect to substrate and base), Figure42. Regioselectivity of reaction with asymmetric substrates follows Zaitsev's rule[206] which states that in the case of possibility of several alkene formation the one with the least number of hydrogens on double bond is formed. Moreover, this reaction is stereoselective: leaving group and proton should be located in antiperiplanar position for effective elimination. The latter explains why different diastereomers form products with opposite orientation of substituents of the double bond.

Figure 42 .

 42 Figure 42. Mechanism of E2 reaction.

  Scheme 4.

Figure 43 .

 43 Figure 43. Example of endo (top) and exo (bottom) cycloaddition. Notice that the products of the reactions are diastereomers.

  External test set was selected for evaluation of the model and comparison with quantum chemistry calculations. As the result 739 equilibria in different conditions were selected for training set and 46 equilibria were selected to external test set. Additional comparison of the model with commercial Tautomerizer tool (ChemAxon) for prediction of tautomer populations in water under room temperature was done. The results of modeling, quantum chemistry benchmarking and commercial tool comparison were published in the article in Journal of Computer-Aided Molecular Design that is shown below. The last modeling challenge was modeling of biorthogonal reactions. The project was initiated within a collaboration with Prof. Frederic Taran, CEA. The initial idea was to develop a model to predict reaction rate for dipolar cycloaddition of sydnones to strained alkynes. The previous Diels-Alder model was useless for this set, as it does not contain reactions of 1,3 dipolar cycloaddition. The training set provided by group of Prof.

Figure 44 .

 44 Figure 44. Studied reactions of sydnones with strained alkynes.

  Moreover, in[233], it was shown that the second step on reaction pathway of sydnone-strained alkyne cycloaddition has almost no barrier which contradicts with kinetic experiments of Liu et al[6]. However, in these two works different sydnones and alkynes were studied. Potentially, fast fluoro-N-phenylsydnone reaction with BCN might have completely different reaction path than medium-rate reaction of unsubstituted Nphenylsydnone. In this work we performed QSPR modeling and quantum chemical study of fluorosydnone reaction pathway.

Fukui

  indices are used to characterize atom ability to share/withdraw electronic charge. For calculation of Fukui nucleophilicity 𝐹 °K, elecrophilicity 𝐹 °± and radical attack susceptibility 𝐹 °² indices of atom A single point calculations of molecules with added and removed electron were done: 𝐹 °K = 𝑃 °(𝑁) -𝑃 °(𝑁 -1) 𝐹 °± = 𝑃 °(𝑁 + 1) -𝑃 °(𝑁) 𝐹 °² = 0.5 * 𝑃 °(𝑁 + 1) -𝑃 °(𝑁 -1) where 𝑃 °(𝑀) -Hirshfield charge on the atom A in molecule with M electrons, 𝑁number of electrons in neutral molecule. Geometry of cation and anion-radical molecules were approximated to the one corresponding to the lowest energy structure of the neutral molecule.

  procedure included scanning along reaction coordinate, optimization of transition state, intrinsic coordinate following to unite transition state with reagents and products. Energetic profile of reaction is shown on Figure 47. N1-O5 length could serve a good approximation for reaction coordinate, it smoothly elongates along the reaction path: 1.379 Å in reagents, 1.414 Å in TS1, 1.546 Å in intermediate and 1.721 Å in TS2 of reaction 15.

Figure 47 .

 47 Figure 47. Reaction pathways for reaction 2 (dashed red line), 15 (solid black line) and 16 (dashed blue line). Relative free energies at 298K of molecules with respect to reagents are shown. Structure of reagents, transition states and intermediate for reaction 15 is shown. Substituent R of sydnone and almost all atoms of BCN are omitted for the sake of clarity. Bond orders correspond to molecule representation on Figure 44.

  3.3.2). The approach includes (i) conformational sampling of initial compounds, (ii) subsequent elimination of local minimum geometries by semi-empiric and DFT calculations, (iii) pre-orientation of reagents to form a valid guess for TS1, (iv) energy screening along the reaction coordinate, (v) optimization of TS1, (vi) selection of lowest energy structure of TS1, (vii) manual examination of TS1 structure, correction and recalculation if required.

Figure 48 .

 48 Figure 48. Free energies of activations predicted by quantum chemical calculations and calculated from experimental values using Eyring formula. Numbers corresponds to reactions in Table 7. Dashed line corresponds to perfect prediction, solid line represent linear correlation.

Figure 49 .

 49 Figure 49. Comparison of HOMO and LUMO energies and symmetries for reactions 2, 15, 16. HOMO and LUMO energies are obtained using DFT scheme on PBE/3z (TZVP) level.

Figure 50 .

 50 Figure 50. Dependency of rate constant on (a) LUMO energy of sydnones and (b) Fukui nucleophilicity on C3 atom. Reaction numbers corresponding toTable 7 is shown next to points.

Figure 51 .

 51 Figure 51. Dependency of rate constant of congeneric series of halo-sydnones on (a) charge of C3 atom of sydnone and (b) LUMO of sydnone. Correlations are shown as dashed lines. Red point and lines correspond to reactions with TMTH, blue ones to reactions with BCN.

Figure 52 .

 52 Figure 52. Transition state of cycloaddition reaction between sydnone 1b (X=I) and BCN. Only few atoms of both reactants are shown for clarity. Dashed line shows steric contacts H…I between halogen atom of sydnone and H atoms attached to one of carbons at the triple bond of alkyne, dotted lines indicate new single bonds forming in this reaction.

3 .

 3 The workflow for reaction property modeling was proposed. It incorporates usage of Condensed Graph of Reaction based fragment descriptors for encoding chemical transformation in combination with solvent and temperature descriptors to represent reaction conditions. Having descriptor vector one can use any machine learning method for model creation (SVM and RF were tried). Using developed approach the model predicting the rate or equilibrium constants of reactions involving various reagents, which occur in many organic solvents and water-organic mixtures were built the first time. It is shown that the RMSE of prediction is comparable with the level of experimental noise.The analysis of prediction errors also shows that the quality of the model is sufficiently high for the identification of data errors and objects with the unique structure with respect to this set of reactions. 4. New type of descriptors for chemical reactions based on mixture representations of reagents and products using SiRMS approach was proposed. Different types of structural descriptors of chemical reactions were benchmarked on cleaned E2 reaction dataset. It showed that for this particular dataset three best structural descriptors are SiRMS mixture descriptors, CGR-based ISIDA fragment descriptors and Morgan fingerprint-based difference reaction fingerprint. 5. The study has clearly shown the importance of correct validation scheme for unbiased estimation of predictive performance of chemical reaction. Two different strategies that are superior to classic cross-validation were proposed: the first one based on calculating predictive performance metrics only on point for which reaction property was measured in only one condition, the second one based on stratified product-out crossvalidation. It was shown that these strategies avoid too optimistic estimation of models performance.

Table 1 .

 1 Données expérimentales utilisées dans ce travail

	réactions	Types de	collectées	Données	standardisées	Données	Jeu modèle	modèle	dans le jeu	Transformation	bibliographiques	Sources
		SN2	7848	7544	4830		1382	[1]
		E2	1431	1389	1043		843		[1]
		CA	1178	1130	880		679		PhD thesis
											defended
											in KFU
	TAU	905	840		782		367		[1]

Table 2 .

 2 Paramètres des jeux de données et performances des modèles a Pour chaque jeu de données, les paramètres suivants sont donnés : le nombre de données (Ntot), le nombre de réactions (Nreact), le nombre de réactions "uniques" (Nunique) pour lesquelles seulement une mesure est disponible, et les coefficients de détermination de validation croisée (R 2 et Runique 2 ) et les erreurs quadratiques moyennes (RMSE and RMSEunique), estimés, respectivement, pour le jeu de données entier et pour son sous-ensembe de données "uniques".

	Dataset	Ntot	Nreact	Nunique RMSE RMSEunique R 2	Runique 2
	SN2	4830	1382	554	0.37	0.65	0.9	0.73
	E2	1043	843	395	0.72	0.87	0.75 0.58
	DA	880	679	279	0.95	1.37	0.8	0.58
	TAU	782	367	267	0.74	0.94	0.79 0.84

a

Table 3 .

 3 Prédiction des constantes d'équilibre tautomérique (logKT) en solution.

	Performance des calculs DFT et des modèles SVR appliqués à deux jeux de
	données externes TEST1 et TEST2. a
	Method	Dataset	Neq	RMSE	R 2	MT,%
	DFT	TEST1 TEST2	23 24	1.62 5.8	0.5 -1.2	74 48
	SVR	TEST1 TEST2	23 24	1.2 2.6	0.73 0.45	74 52

  Quantum chemistry is still the only one way for ab initio prediction (e.g. not based on previous knowledge of related examples) of energetic characteristics like rate constants or activation barriers. But it is limited by intrinsic inaccuracy of polarizable continuum solvent models. Additionally, fast methods like DFT must first be calibrated on certain sets of molecules and usually does not achieve "chemical" accuracy (about 1 kcal/mol)[START_REF] Peverati | Quest for a universal density functional: the accuracy of density functionals across a broad spectrum of databases in chemistry and physics[END_REF]. The accuracy versus computing effort ratio is much too low to qualify these approaches as large-scale predictors of properties of large series of reactions. So, we decided to use another strategy, which is empirical, not mechanistic and, in exchange, extremely fast. The proposed approach implies application of QSAR/QSPR modeling, a

The main advantage of this method is that the whole energetic profile of a reaction can be generated by calculations. It means that in principle all possible ways of reaction can be determined with rather high detail resolution of the electron transferring processes. machine-learning technique that learns a model of a property from known examples and extrapolates property values of new items based on this model. The Condensed Graph of Reaction (CGR) approach helped exploiting rich base of Chemoinformatics methods to predict properties of reactions. CGR representation (as it shown at

  networks. Dataset comprised 3451 SN2 reaction rate constant in pure solvents at different temperature. Authors used fragment, topological and quantum-chemical descriptors for encoding chemical transformation, Palm solvent descriptors to represent media. descriptors that encode concentration of substrate and nucleophile and the nature of counter ion increases predictive performance of the model.In this case RMSE achieved 0.98 vs 1.07 logk units in case of absence of counter ion descriptor. Such wise despite 4 models for SN2 reactions were published by now, there was no model that was developed to predict reaction rate of any nucleophile in any solvent and solvent mixtures.Thus in the work ofSchneider et al [100] several machine learning techniques were used to make classification model trained on dataset of reactions, extracted from patents and annotated by type using NameRxn tool. Information on conditions could help with classification since reactions of a given type are usually performed in similar conditions.

	Thus, authors used very unusual fingerprints: they substracted from product fingerprint
	reagent one and summed result with fingerprint of small-molecule reactant. Several
	2.4.2.2. RDKit molecular fingerprints were used for reaction fingerprint creation [148]. The Reaction classification
	Classification of reaction types implies the recognition of types of reactions present created model were used to annotate reactions that were skipped by rule-based approach.
	in databases or entered by a user, which is important for solving practical tasks of
	synthetic chemistry, for example, for search of similar reaction, optimal condition
	Data-based classification approaches are based on application of machine learning
	tools to the dataset with known reaction classes to create the model able to classify new
	was built using ISIDA and SiRMS descriptors of structure, 15 solvent descriptors reactions. One of first approaches to classify reaction by type have used Kohonen's map
	characterizing polarity, polarizability, and proton-acceptor and proton-donor abilities of and physicochemical descriptors of atoms [146]. It was shown that reactions could be
	solvent, temperature of reaction. Random Forest regressor usage for model building quite efficiently classified; however the number of reaction types was quite small.
	allowed reaching RMSE on out-of-bag sample of 0.5 logk units that was shown to be Similarity based system for hierarchical reaction classification to additions, eliminations,

Significant descriptors were selected using the fast stepwise linear regression procedure.

On 193 test set reactions RMSE was 0.58 log units. The approach was used as well for SN2/ SN1 reaction classification and SN1 reaction rate prediction

[START_REF] Kravtsov | Prediction of the preferable mechanism of nucleophilic substitution at saturated carbon atom and prognosis of S N 1 rate constants by means of QSPR[END_REF]

with RMSE 0.61 log units on test set that is unbiased enough.

Further development of reaction rate modeling was related with CGR approach.

Hoonakker et al [139]

used dataset of 1014 reactions in water solvent proceeding at various temperatures, 3 various machine-learning techniques (SVR, M5P, MLR) and ISIDA fragment descriptors based on CGR representation of reactions for building model for logk prediction. The best model built using SVR and atomic sequences of topological length from 2 to 8 as descriptors shown quite good performance, Q 2 =0.53, RMSE=1.26, but worse than that of Kravtsov et al

[START_REF] Kravtsov | Prediction of Rate Constants of SN2 Reactions by the Multicomponent QSPR Method[END_REF]

. The difference could arise due to different validation procedures and biased estimation of prediction error (see chapter 4.2). Hoonakker's model could only predict reaction rate in water solvent. In the next work

[140]

the dataset of 1041 reaction proceeding in different solvents and at different temperatures involving only neutral nucleophiles was collected. The model on this data comparable with experimental noise. Neutral nucleophiles were selected to exclude effects due to complication of reaction process in case of usage of salts as nucleophiles. In the next work

[141]

it was found that for particular case of SN2 reaction using azide ion as nucleophile addition of selection or synthesis planning (generation or selection of retrosynthetic rules). There is two approaches for reaction classification: model-based and data-based

[START_REF] Chen | Reaction Classification and Knowledge Acquisition[END_REF]

. Modelbased approaches use different reaction center representations schemes (see Chapter 2.2.2). Reaction center is specific for reaction of given type and that could be utilized for grouping, making reaction nomenclature and ontology

[START_REF] Hendrickson | The Variety of Thermal Pericyclic Reactions[END_REF][START_REF] Hendrickson | Comprehensive System for Classification and Nomenclature of Organic Reactions[END_REF]

. It is however not a universal approach since classification of reactions used in organic chemistry sometimes is not related to immediate environment of reactions center.

Similar approach that could be called rule-based works in opposite way: first, reaction ontology is manually created using concepts in organic chemistry and then rules for reaction extraction are manually proposed. For example, RXNO reaction ontology (http://github.com/rsc-ontologies/rxno) created on the basis of works

[142, 143]

was used in rule-based reaction classification NameRxn system [144]. This system was used to investigate the popularity of reactions of various types in medicinal chemistry

[145]

. and substitutions, followed by two successive subdivisions by number and types of reactive atoms was proposed by

Sello and Termini [147]

. Rules in mentioned NameRxn system lack generality and thus they sometimes skip reactions belonging to given type.

  The model was built using manually annotated set of 193 reactions.

	Several descriptor types (reagent-based, difference fingerprints, CGR-based fragments)
	and machine learning techniques were benchmarked. The best model had balanced
	accuracy 0.85± 0.15. The model was successfully tested on external set and was capable
	to predict feasibility conditions for new 50 Michael reactions.

  Chemical reactions study and discovery is a central topic in synthetic and quantum chemistry. But as it comes from the review reaction modeling is really poorly explored topic in chemoinformatics. There are few sporadic publications without systematic efforts for the development of approach adapted for reaction modeling taking into account reaction specifics. The goals of present study immediately come from previous studies. In

	this PhD thesis we want to contribute to reaction modeling by (i) collection of big dataset
	of kinetic and thermodynamic properties of chemical reactions, (ii) development of
	workflow for modeling reaction characteristics using machine learning techniques, (iii)
	development of specific methodologies for reaction data curation, (iv) development of
	novel descriptors to model chemical reactions taking into account reaction conditions, (v)

The authors predicted yields and times of reactions with several types of descriptors. The author claimed that the results of their work were somewhat negative but tended to be thought-provoking. Error in prediction of yields and reaction times are 35% and 25%, correspondingly. Classification models have accuracy 65% for prediction of reactions of class with high yields (>65%)

[150]

. Such limited prediction was interpreted as consequence of imperfectness of descriptors but one should take into account that yield of reaction is a very noisy parameter.

In the work of

Ahneman et al [151]

high-throughput reaction screening approach for collection data on yield of Buchwald-Hartwig reaction was presented. The yield was predicted using different machine-learning techniques and quantum-chemical descriptors. Unlike aforementioned work, model shows very good performance on out-ofbag validation: RMSE was substantially lower than in previous work -11% and R 2 =0.91. Such a drastic difference between two articles could be partially explained by different dataset (big and noisy data from Reaxys vs relatively small automatically collected data). modeling reactions studied previously and a new one , creation of universal models predicting reaction characteristics of different types in variety of conditions, (vi) development of server for publication of reaction models.

  NumPy (basic math and statistics)[178], ChemAxon Web Services (chemical data processing)[179]. Client side for model publishing server was written in JavaScript.Using aforementioned tools several new tools and libraries were created. For databasing chemical reactions in-house chemical cartridge CGR DB was implemented.

	CGRtools library were created for reaction information management using CGR
	approach. For QSRR modeling CIMMtools library was implemented. All the
	developments were supervised by Dr. R. Nugmanov and Dr. T. Madzhidov from
	Chemoinformatics and Molecular Modelling Laboratory of KFU. The final code was
	implemented by Dr. R. Nugmanov on the base techniques that were developed during
	this work.
	5 language [172]. For
	handling graphs NetworkX library [173] is used. Sklearn [174] library is used for
	machine learning method (mainly SVM) application. Databases are implemented using

Pony ORM

[175]

solution linked to PostgreSQL database management system [176]. Other dependencies include the following libraries as well: Pandas (data operation) [177], GTM method for visualization and modeling as well as the Fragmentor program for fragment descriptor calculation and Genetic Algorithm optimizer for model hyperparameter selection are the developments of the Laboratory of Chemoinformatics of the University of Strasbourg, supervised by Prof. A. Varnek.

  The CGR looks like pseudomolecule, this similarity was used to store CGRs in SDF file like an ordinary molecule. The solution is compatible It is shown from the first field (starting with "M SAL") which atom or bond CGR property corresponds to. Precisely first field encode atomic array that can contain 1 entry (for dynamic atom) or 2 entries (for dynamic bonds). The second line

	with a standard MDL/SDF format. Common bonds are encoded according to
	MDL/SDF specification. Dynamic bonds are encoded additionally in bond section as 8
	("any bond" according to the CTfile specification [181]).
	The type of dynamic property label is encoded in MOL file properties section using

standard V2000 connection table specification of specific groups (Sgroups) in CTfile format

[181]

, see Scheme 1 below. First line in properties section starting "M STY" in property section specifies number of dynamic properties. Then, each dynamic property is encoded in 4 fields.

  Transformation table contain structural representation of reaction transformation, its CGRS and CSMIRKS signatures. As additional fields MD5 hashed representation [194] of CSMIRKS and CGRS were used. The transformation table has one-to-many relation to conditions table. Hashed CGRS and CSMIRKS representation are used for fast reaction search and finding transformation for a given reaction. The latter is needed to establish relations in database.

	•	Conditions table contain information about conditions used for
	measurement and corresponding reaction property.

Table 5 .

 5 Data that were collected and curated for this research.

	Dataset	collected	Data	Curated data	Model set	set	ons in model	Transformati	data	Source of
	SN2	7848	7544	4830	1382		[1]
	E2	1431	1389	1043	843		[1]
	CA	1178	1130	880	679		PhD thesis
									defended in
									KFU
	TAU	905	840	782	367		[1]
	Every reaction/equilibria has four layers of associated information:

•

Transformation -description of structural changes in reaction/ equilibria: o structures of reagents and products in MDL RDF format(two forms of tautomers were stored as reagent and product), o atom-to-atom mapping, based on reaction type reported in the source.

• Conditions -description of media and physical conditions at which reaction rate constant was measured o temperature (in Celsius)

o solvent name (only organic solvent for solvent mixtures with water)

  , see Chapter 2.4.3.1. However none of them could

	predict reactions in different solvents including water-organic solvent mixture and with
	different types of reagents.
	NH 2	H 3 C	NH 2 +
	Br CH 3		Br
	H 3 C	H 3 C

  This is a case of reactions 14 and 15, which differ only in solvent. Their rate constants differ by two powers of magnitude. When one of these reactions is selected to test set, the model predicts logk shifted towards reaction rate of the training set reaction.

		Table 6. Experimental ("exp") and predicted ("pred") rate constant logarithms for
	nucleophilic substitution reactions involving anionic nucleophiles.
	N					Reaction						Conditions Exp Pred
								CH 3
	Cl
			+	H 3 C	OH N	CH 3		H 3 C	N +		OH	+	Cl -	Phenyl-ethanol 100 %, 201 °C	-3.21 -1.36
	O	N +	O -
									O	N +	O -
		I	Br	H 3 C	O -	Na +	Br	O CH 3		I -		Na +	methanol 100 %, 50 °C	-4.99 -3.08
												methanol 100 %, 0 °C	-1.58 -3.47
												toluene 100 %, 30 °C	-3.30 -5.98
												methanol 100 %,0 °C	-5.00 -2.85
												methanol 100 %,55 °C	-5.70	-3.46
												DMFA 100 %, -20 °C	-0.33 -3.08
												water 100 %, 25 °C	-5.80 -2.38
												toluene 100 %, 30 °C	-2.65 -4.53

  So, local models behave like global one and splitting of the dataset in two separate does not provide any rise in accuracy.

	I CH 3 model)=0.59.	+		OH -	+		Na +			OH	+	I -	+	Na +	DMSO 100 %, 25 °C	-2.48	1.33
		HO	HO	Cl	+		OH -	+	Na +		HO	HO			OH	+	Cl -	+	Na +	water 100 %, 0 °C	-2.06 -6.38
	I	CH 3	+	O	O -CH 3	+	H 3 C	H 3 C	N +		CH 3 CH 3	H 3 C O	O CH 3		+	I -	+	H 3 C	H 3 C	N +	CH 3 CH 3	DMFA 100 %, 0 °C	0.29	-1.90
														N				N +
								I CH 3	+											+	I -	toluene 100 %, 30 °C	-4.01 -6.04
													-O	N +	O	O		N +	O -
	I	CH 3	+	O	CH 3 O -	+					As +			H 3 C O	CH 3 O	+	I -	+	As +	acetonitrile 100 %, 25 °C	-0.32 -2.71
	I	CH 3	+	O	CH 3 O -	+					As +			H 3 C O	CH 3 O	+	I -	+	As +	metanol 100 %, 25 °C	-5.57 -3.26
						Cl															OH
					O													O
					S	O		+		OH -	+	Na +					S	O		+	Cl -	+	Na +	1,4-dioxane 100 %, 50 °C	-2.05 -4.40
		Br		CH 3 +	-O	O N + O -		+	Li +		H 3 C		O		O O -N +	+	Br -	+	Li +	ethanol 100 %, 25 °C	-7.48 -5.31
					I		I	+	H 3 C	O -	+	Na +	I	O CH 3	+	I -	+	Na +	methanol 100 %, 20 °C	-7.00 -4.73
		5.6.		Local and global models
		Results reported previous section reveal that model's performance assessed on the
	reactions with neutral and anionic nucleophiles is similar. Thus, the following parameters
	were obtained in cross-validation: R 2 =0.91 and 0.91, RMSE=0.39 and 0.3 for anionic
	and neutral nucleophiles, respectively. Similar situation is observed for a subset of unique
	datapoints: R 2 UDP = 0.75 and 0.74, RMSEUDP = 0.68 and 0.57 for anionic and neutral
	nucleophiles, respectively. So, it was interesting to build models for datasets, containing
	Br -only particular types of reactions (here, with neutral or anionic nucleophiles), we call Ag + + + + +
	N them local models to distinguish from global models built on the entire set. According to H 3 C CH 3 CH 3 CH 3 N + H 3 C UDP-based validation local model performance is very close to the global model one
	CH 3 shown above, RMSEUDP (anion, local model) = 0.72 and RMSEUDP(neutral, local I I -+ +
														O	CH 3			O	CH 3

  The ability to selectively form and break chemical bonds in chemically complex and uncontrollable biological media is a long-standing goal of chemists interested in modifying biological materials. Biorthogonal chemical reactions[222], which are reactions that do not interfere with biological processes, address precisely this challenge and therefore they are of major importance in the fields of chemical biology and biochemistry. To fulfill the requirements of bioorthogonality, reaction partners must be stable and inert towards the plethora of chemical functionalities found in living systems while reacting selectively, efficiently and rapidly with each other under physiological conditions with no or innocuous by-products.

	Biorthogonal reactions play key roles in modern biochemistry [223, 224]. 1,3-
	dipolar cycloaddition reactions like copper-catalyzed azide-alkyne cycloaddition
	(CuAAC) are one of the most extensively studied and used biorthogonal transformations.
	But they were catalyzed by Cu salts, which resulted in a dramatic toxicity of this systems
	for living organisms [225]. New reaction systems based on strain promoted cycloadditions
	were explored [226]. The catalyst addition is not required for them but reaction have
	rather low rate [227].
	Recently sydnones were shown to be a very promising class of biorthogonal 1,3-
	diene that can participate cyclization reaction with cyclic alkynes resulting in pyrazole
	formation [228]. This reaction proceeds in two step: cyclization with formation of
	reactive intermediate and retro-Diels-Alder reaction with CO2 cleavage (Figure 44).

Thermal cycloaddition of sydnones with alkynes require quite harsh conditions and proceeds with low regioselectivity

[228]

while under copper catalysis reaction proceed in mild conditions with high yield and selectivity

[229][230][231]

. Strained alkynes readily react with sydnones and iminosydnones without catalyst

[232][233][234]

. Substituents in sydnone ring

Table 7 .

 7 Dataset of sydnone-alkyne strain-promoted cycloaddition reaction rate constants. R-group position are shown on Figure 44.

	No.	R	X	Cycloalkyne Rate	constant	Reference
					(tolerance), M -1 *sec -
					1
	1	p-MeO	H	BCN	0.006 (0.001)		[235]
		C6H4
	2	p-Me C6H4 H	BCN	0.032 (0.001)		[6]
	3	C6H5	H	BCN	0.027 (0.002)		[235]
	4	p-CO2H	H	BCN	0.059 (0.001)		[235]
		C6H4
	5	p-CF3 C6H4 H	BCN	0.199 (0.002)		[235]
	6	p-NO2	H	BCN	0.289 (0.012)		[235]
		C6H4
	7	C6H5	CH3	BCN	0.018 (0.002)		[235]
	8	C6H5	C6H5	BCN	0.027 (0.001)		[235]
	9	C6H5	CF3	BCN	0.008 (0.001)		[235]
	10	p-Me C6H4 Cl	BCN	0.872 (0.034)		[6]
	11	p-Me C6H4 Br	BCN	0.592 (0.021)		[6]
	12	p-Me C6H4 I	BCN	0.306 (0.008)		[6]
	13	p-CO2H	Br	BCN	0.798 (0.065)		[235]
		C6H4
	14	p-CO2H	Cl	BCN	1.593 (0.034)		[235]
		C6H4
	15	p-Me C6H4 F	BCN	42		[6]
	16	p-Me C6H4 F	TMTH	1500		[6]
	17	p-F C6H4	F	TMTH	3500		[6]
	18	p-CF3 C6H4 F	TMTH	12000		[6]
		From the table one can see that reactions with TMTH have several magnitudes
	higher speed than with BCN and their rate were measured only for the most reactive
	fluorosydnone. We removed from dataset 3 reactions, that have the same reaction rates in

two articles, but have difference in structure (reactions

10,11,12 in [234]

have CH3 group in para position of phenyl group, but in previous article

[229]

they have not). So, out of six reactions in the article

[229]

we kept only 3.

3.3. Activation free energy calculation

  Pair of enantiomeric structures of transition states should have the same energies.However, due to some fluctuations, structures of them obtained using 4 different orientations of reagents usually had slightly different energies. For further analysis the lowest energy transition state among four was selected. The selected structure hessian was calculated and the presence of a single large imaginary frequency was checked. The correctness of transition state was checked using Intrinsic Reactional Coordinate following (IRC) procedure.

	Free energy of molecule and transition state formation was calculated using built-in
	thermochemical calculation in rigid ideal rotator and harmonic oscillator approximation
	at 298 K in Priroda11 program, including zero-point vibrational energy corrections
	[165].
	Formula Eq. 6 was used for calculation of activation free energy on the basis of
	quantum chemical calculations:
	Eq. 6	∆∆G g = ∆G =›oe `• -(∆G =›oe žŸ7C6C@ +∆G =›oe A ^ŸC@ )
	Calculation of activation free energy based on experimentally measured reaction
	rate was done using formula Eq. 7 based on Eyring equation of Transition state theory
	Eq. 8:
	Eq. 7	𝑘 = 𝜅	¢ _	à 𝑒	£∆∆¤ e NO
	Eq. 8	∆∆𝐺 g = -𝑅𝑇ln (	^ * a ¥¢ _ `)
	Relativistic effects were neglected since for molecules under study it has minor
	importance. Priroda11 is probably one of the fastest DFT code due to efficient evaluation
	of density functional exchange-correlation terms based on the expansion of the electron
	density over auxiliary basis set [166]. The program was used since computational
	efficiency was of major importance for the study.
	Geometry optimization for reagents, intermediates and products as well as saddle
	point optimization were performed using built-in quasi-Newton-Raphson procedure and
	BFGS hessian update scheme. Scanning along plausible reaction coordinate was used to
	localize good structural guess for transition state. The scanning procedure was
	constrained local optimization with one coordinate N1…Calkyne set externally. Its value
	changed from 2.1 to 3.1 with step size 0.1 Å. Geometry optimization was followed by
	frequency calculation to control correct structure of hessian: discussed structures of
	reagents, products and intermediate had no imaginary frequency, transition state

calculations were made in

Priroda 11 program [165]

with PBE exchange and correlation functional

[126]

and built-in triplezeta split valence basis set (called 3z, equivalent of

Schäfer's TZVP basis [167]

). Figure 46. Four possible pre-reaction complexes diene and dienophile 9.

where k -rate constant; κ -transmission coefficient (here, κ=1); T -Temperature in Kelvin; KB -Boltzmann constant; ∆∆G ≠ -Difference of free energies of reagents and TS at 298.15 K; h -Planck constant; R -gas constant.

9.3.4. Conceptual DFT indices

Different Conceptual

DFT [244]

indices were used in the work to reveal structural factors responsible for reaction rate. Electrophilicity index, ω, measures the stabilization energy when system acquires an additional electronic charge [246]:

Acknowledgements I would like to express my gratitude to my academic advisors, Prof Alexandre Varnek and Dr Timur Madzhidov for their help in academic challenges that I have faced. I would like to thank all the members of the jury, Dr Michel Petitjean, Dr Peter Ertl and Prof Esther Kellenberger for accepting to judge and review my work. I would also like to express my sincere gratitude to all colleagues and co-authors from both laboratories of Chemoinformatics from Kazan and Strasbourg: Dr. Dragos Horvath, Dr. Gilles Marcou, Dr. Pavel Polishchuk, Dr. Ramil Nugmanov, Dr. Fanny Bonachera as well as, former and current PhD students Dr. Helena Gaspar, Grace Delouis, Pavel Sidorov, Marta Glavatskih, Arkadii Lin, Iuri Casciuc for their help in the work, advice and support. Special thanks to Dr. Olga Klimchuk and Soumia Hnini for their help on accommodation and administration tasks. I am thankful to all people that helped to collect experimental data from the literature in Kazan Federal University: N. Khafizov, N. Shalin, A. Bodrov, D. Malakhova, A. Petrovski, F. Bekmuratova, G. Khaliullina, and some others.

Chapter 3.

Computational techniques used in the study

implemented by Laikov [166]. PBE exchange-correlation functional was used [126].

Built-in triple-zeta split valence basis set (called 3z within program itself, equivalent of Schäfer's TZVP basis [167]) was selected. This level of theory was used for geometry optimization, transition state localization, reaction path exploration and Intrisic Reaction Coordinate following along. Transition states were found with saddle point search algorithm of Priroda program.

Gaussian program [168] was used for solvation free energy calculation. In this case IEF-PCM model [START_REF] Tomasi | The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level[END_REF]169] with SMD parameters for non-electrostatic terms [170] was utilized. Geometry optimization and Hessian calculation in solvent was performed in Gaussian program [171] using PBEPBE/6-311++G(d,p) level.

All structures under discussion were optimized using default Gaussian and Priroda program algorithms. Many different starting geometries were used in order to enhance chance to find global minima of energy. Frequency calculation supported that discussed structures have the right set of Hessian eigenvalues: all positive frequencies for molecules and intermediates, one imaginary frequency for transition states. The details are shown in the chapter devoted to 4.3.3.

3.4.

Model implementation

For model publication, a client-server application was implemented. The description of interaction with client side is described in this chapter.

All QSRR models are available at cimm.kpfu.ru/predictor after registration procedure. The first page of predictor contains User manual and tools for interactive creation of molecules or reactions (Add task, Figure 23) and button for files upload. After creation of some dataset one should click the button "Validate" to check structure for correctness.

Chapter 4.

Data collection, cleaning and representation

Chapter 5.

Models for rate constants of bimolecular nucleophilic substitution reactions

Chapter 6.

Modeling of rate constants of bimolecular elimination (E2) reactions modeling

Conclusive remarks

The first model predicting the rate constant of E2 reaction with different reagents and in wide variety of conditions was prepared using fragment descriptors generated for Condensed Graph of Reaction and some special descriptors accounting for experimental conditions. The predictive performance of the consensus model in cross-validation was high enough: RMSE = 0.69 log units, R 2 = 0.75. Analysis of outlier was shown that model fails to predict logk for reactions involving substrates with unique fragments . The model was published online. The results supported the universality of the approach.

Models built on mixture-based reaction descriptors

The question rises how well the CGR-based descriptors perform in comparison with other known descriptor types for reaction. Thus benchmarking study of different descriptor types was made. CGR-based fragment descriptors were compared with several difference fingerprints calculated by RDKit and proposed in the work [100]. Moreover the novel type of descriptors suitable for reaction modeling was introduced. They are based on the SiRMS descriptors for compounds mixture [START_REF] Oprisiu | QSPR Approach to Predict Nonadditive Properties of Mixtures. Application to Bubble Point Temperatures of Binary Mixtures of Liquids[END_REF]. Formally left-and righthand side of reaction could be considered as mixture of reactants and products respectively and corresponding SiRMS mixture descriptors for them could be calculated.

Resulting descriptor vector could be combined by concatenation or subtraction. Hence three completely orthogonal descriptor generation strategies were compared: difference fingerprints, CGR-based fragment descriptors and mixture descriptors for chemical reaction.

Random Forest [153] (RF) machine learning algorithm for model building. For the study we needed non-linear regressor efficiently working with very large descriptor space (when number of descriptors is substantially larger than number of reactions to model), having as few hyperparameters as possible to adjust. RF perfectly suited our need due to efficient tackling non-linearities, non-sensitivity to descriptor vector size, and only one hyperparameter that is required to be adjusted -ratio of descriptors that is randomly selected for tree branching (number of trees in the forest should be as large as possible, we used 500, and the other parameters of RF influence poorly and default values could be accepted).

Conclusive remarks

In this section, we described new type of descriptors for chemical reaction based on mixture representation of reagents and products. Unlike CGR-based descriptors this approach does not require Atom-to -Atom Mapping as a data pre-processing step. Reaction-out cross-validation and product-out cross-validation lead to quite different ranking of descriptors type according to the models performances. The bias introduced in reaction-out cross-validation could drastically affect the conclusions drawn on models performance.

Chapter 7.

Modeling of rate constants of Diels-Alder reactions

Conclusions remarks

The CGR-based approach and developed workflow were used to build consensus model for Diels-Alder reactions rate constant prediction. The model displays a reasonable predictive performance: R 2 =0.87 and RMSE = 0.75 log units in cross-validation. The model was published on line on http://cimm.kpfu.ru/predictor.

Conclusive remarks

In this section we reported the first attempt to model tautomerisation equilibrium Along with global model built on whole training set, local models built on particular type of tautomerisation were tested. In most of cases, local model lead to lower RMSE than the global model. However some tautomerism types are predicted much worse than others. This can be explained by small number of datapoints and by its structural diversity.

The comparison with other methods as ChemAxon Tautomerizer plugin or quantum chemistry calculations showed that developed approach has higher precision in prediction on external test.

The model was published on-line on the site http://cimm.kpfu.ru/predictor.

Summing up the results shown in the Chapters 5-8, we could conclude that the efficiency of CGR-based and solvent descriptors in modeling of reaction characteristics.

Prediction errors were found at the level of the data noise. The developed approach was used to build models to predict rates of reaction of different classes: substitution (SN2), elimination (E2), cycloaddition (Diels-Alder) reactions, which were built for the first time;

as well as tautomeric equilibrium constant model utilizing new methodology of direct prediction without intermediate building of acidity models.

Chapter 9.

Modeling of reaction rates of some bioorthogonal reactions

geometry had one large imaginary frequency. The correctness of transition states was also checked by intrinsic reaction coordinate (IRC) following approach.

For 3 reactions full reaction paths were calculated. First stage transition state (TS1) was localized and optimized using aforementioned procedure. in the other direction it converges to the final product molecules.

Solvation free energy was calculated for some structures by IEF-PCM model [START_REF] Tomasi | The IEF version of the PCM solvation method: an overview of a new method addressed to study molecular solutes at the QM ab initio level[END_REF]169] with SMD parameters for non-electrostatic terms [170] using Gaussian program [171]. In this case geometry optimization and hessian calculation was performed in

Gaussian program [171] too using PBEPBE/6-311++G(d,p) method.

Procedure for transition state detection

The workflow used for semi-automatic first transition state detection is schematically represented on Figure 45. The detailed descriptions of each step are given below.

Initial conformers generation and optimization using MMFF94 force field

Initial conformers generation

All structures were drawn following the same numeration of sydnone ring atoms.

ChemAxon Calculator Plugin cxcalc was used to generate up to 500 conformers with diversity limit 0.1 Å with subsequent optimization using MMFF94 force field.

PM6 optimization

All structures obtained from previous step were optimized using PM6 semi-empiric method using MOPAC 11 program [245]. Duplicated conformers (with RMSD <0.1 Å) were removed using simple Perl script that aligns molecules according to principal components of inertia with following RMSD calculation (developed at N.N.Vorozhtsov

Institute

Reagent structure optimization using DFT

All structures obtained from PM6 optimization were optimized in PBE/3z in Priroda11. The lowest energy conformer was selected to continue with the next step. The selected conformer geometry was optimized until no frequency of the hessian was imaginary.

Pre-reaction complexes formation

There are 4 possible orientations of reagents to form TS1, that could be considered as 2 pairs of enantiomers, Figure 46. However when CO2 is cleaved the product is formed as racemic mixture of two enantiomers.

Four possible orientations of reagents were generated using in-house Python script in a way that distance between reaction center atoms (C3…C1 and N1…C2 or C2…C1 and N1…C2, depending on orientation, atomic numbering is given on Figure 46) was 3.1 Å, i.e. much greater than in transition state. Guess for first transition state was found using scanning procedure in Priroda11 varying distance between reaction center atoms from 3.1 Å to 2.1 Å with step size 0.1 Å. For the structure with the lowest energy hessian was calculated, and saddle point optimization started.

As one can see, quantum chemical calculations in gas phase well describe activation free energy of sydnones cycloaddition to strained alkynes. Root mean squared deviation between predicted and experimental values is 1.97 kcal/mol, the plot predicted vs experimental is shown on Figure 48. Considering that the accuracy of the free energy measures is about 1 kcal/mol, the agreement between the models and the experiments is reasonable and is comparable to what is expected using common DFT functional [START_REF] Peverati | Quest for a universal density functional: the accuracy of density functionals across a broad spectrum of databases in chemistry and physics[END_REF].

Moreover, the calculated value could contain error due to approximate nature of

Structural factors responsible for reaction rate

Quantum chemical assessment of reaction rate based on the developed semiautomatic approach is reliable but computational resource-consuming task since the whole calculation could take some days per CPU core. The most complicated is the localization and optimization of transition state. For real application one needs to reduce search space and find factors responsible for high reactivity. Our attempt to build QSPR models failed since the resulting models have mediocre predictive ability according to cross-validation procedure.

General conclusions