Marta GLAVATSKIKH

Modeling and visualization of complex chemical data using local descriptors

Résumé Cette étude considère des systèmes où non seulement la structure moléculaire, mais les conditions expérimentales sont impliquées. Les structures chimiques ont été codées par des descripteurs locaux ISIDA MA ou ISIDA CGR, ciblant spécifiquement les centres actifs et leur environnement le plus proche. Les descripteurs locaux ont été combinés avec les paramètres spécifiques des conditions expérimentales, codant ainsi un objet chimique particulier. La méthodologie a été appliquée avec succès pour la modélisation QSPR des paramètres thermodynamiques et cinétiques des interactions intermoléculaires (liaisons halogène et hydrogène), des équilibres tautomères et des réactions chimiques (cycloaddition et SN1). La méthode GTM a été appliquée pour la première fois pour la modélisation et la visualisation de données chimiques mixtes.La méthode sépare avec succès les groupes de données à la fois en raison des structures et des conditions.

Résumé en anglais

This work describes original approaches for predictive chemoinformatics modeling of molecular interactions and reactions as a function of the structures of interacting partners and of the chemical environment (experimental conditions). Chemical structures have been encoded by local ISIDA MA-based or CGR-based descriptors, specifically targeting the active centers and their closest environment. The local descriptors have been combined with the specific parameters of experimental conditions, thereby encoding a particular chemical object. The methodology has been successfully applied for QSPR modeling of thermodynamic and kinetic parameters of intermolecular interactions (halogen and hydrogen bonds), tautomeric equilibria and chemical reactions (cycloaddition and SN1). GTM method has been applied for the first time for QSPR modeling and visualization of mixed chemical data. This method successfully separates data clusters on account of both chemical structures and experimental conditions.

Résumé en français

La complexité des données chimiques reste un défi pour la modélisation structurepropriété. En particulier, ceci concerne le développement de modèles prédictifs pour les propriétés liées à des centres sélectionnés (atomes ou groupes chimiques), par exemple les différents types d'interactions intermoléculaires ou de réactions chimiques. Un autre niveau de complexité vient du fait que ces propriétés dépendent souvent non seulement de la structure chimique des molécules en interaction, mais aussi des conditions expérimentales (solvant et température).

Afin de décrire correctement cette complexité, les modèles de propriétéstructure associés doivent impliquer des descripteurs caractérisant à la fois les aspects structurels et conditionnels. De plus, les structures chimiques doivent de préférence être codées par un descripteur local spécial [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] ciblant spécifiquement les centres sélectionnés sur les espèces d'interaction et leur environnement le plus proche.

Cette étude considère des systèmes de niveau de complexité différents dans la plupart desquels non seulement la structure moléculaire, mais aussi les conditions expérimentales jouent un rôle significatif (Tableau 1). Le premier exemple est le plus simple : il concerne la modélisation de la stabilité de la liaison halogène mesurée par la constante d'équilibre de complexes de molécules organiques avec le même accepteur (I 2 ) dans un solvant (hexane) à 298K. Dans ce cas, les conditions expérimentales sont contrôlées et fixes ; seule l'espèce chimique varie.

La complexité augmente avec la modélisation de l'énergie libre de liaisons i. Glavatskikh, M., Madzhidov, T., Solov'ev, V., Marcou, G., Horvath, D., Graton, J., ... & Varnek, A. (2016).

Predictive Models for Halogen-bond Basicity of Binding Sites of Polyfunctional Molecules. Molecular informatics, 35(2), 70-80. ii.

Glavatskikh, M., Madzhidov, T., Solov'ev, V., Marcou, G., Horvath, D., & Varnek, A. (2016). Predictive models for the free energy of hydrogen bonded complexes with single and cooperative hydrogen bonds. Molecular informatics, 35 (11-12), 629-638.

Les modèles ont été construits à l'aide des méthodes SVM (Support Vector Machine), MLR (Multiple Linear Regression) et GTM (Generative Topographic Mapping). La SVM et la MLR sont des procédés d'apprentissage automatique conventionnels largement utilisés en chémoinformatique, tandis que le second, initialement développé comme outil de visualisation de données [START_REF] Bishop | GTM: The generative topographic mapping[END_REF] , a été étendu en laboratoire à des tâches de modélisation structure-propriété [START_REF] Gaspar | GTM-Based QSAR Models and Their Applicability Domains[END_REF][START_REF] Gaspar | Chemical data visualization and analysis with incremental generative topographic mapping: big data challenge[END_REF] . Différents types de descripteurs locaux (à base de MA et à base de CGR, voir section 2) ont été utilisés dans la modélisation de la stabilité des liaisons halogènes, des liaisons hydrogènes et des équilibres de réactions chimiques. Pour ces derniers, les descripteurs structuraux ont été complétés par des descripteurs spécifiques de solvant et de température.

La Quatre scénarios de descripteurs à base de MA ont été considérés 2 :

• Séquences d'atomes à partir de l'atome marqué, avec des fragments centrés sur l'atome avec l'atome marqué central (MA1).

• Fragments contenant des atomes marqués (MA2).

• Fragments avec et sans atomes marqués (MA3).

• N'utilisant pas les atomes marqués (MA0), également utilisés à des fins de comparaison.

La longueur du vecteur des descripteurs varie en fonction de la stratégie sélectionnée. Ainsi, MA0 et MA2 sont des sous-ensembles de MA3, alors que M1 est un sous-ensemble de MA2.

Dans le graphe condensé de réaction (CGR) [START_REF] Hoonakker | Condensed Graph of Reaction: considering a chemical reaction as one single pseudo molecule[END_REF] , les structures de tous les réactifs et produits sont fusionnées en un seul graphe (Figure 1, en bas) qui décrit à la fois les liaisons chimiques conventionnelles (simples, doubles, aromatiques, ...) et dynamiques des liaisons caractérisant des transformations chimiques (par exemple, simple à double, simple brisée, simple créée, etc.). Les fragments contenant des liaisons dynamiques ont été utilisés comme descripteurs locaux.

Les conditions expérimentales ont été codées par 13 paramètres de solvant [START_REF] Madzhidov | Structure-reactivity relationships in terms of the condensed graphs of reactions[END_REF] reflétant la polarité, la polarisabilité, l'acidité et la basicité, et l'inverse de la température

(1/T).

Le vecteur des descripteurs entier pour un processus donné résulte de la concaténation de descripteurs de conditions structurelles et expérimentales.

2. Modélisation quantitative de la relation structure / propriété / réactivité (QSPR) de différents objets chimiques.

Le workflow général de modélisation incluait les étapes suivantes: (1) collecte et conservation de données, (2) génération de différents ensembles de descripteurs ISIDA, (3) sélection du meilleur ensemble de descripteurs en fonction des performances des modèles en validation croisée, (4) construction de modèles individuels et consensus (pour SVM seulement) impliquant des descripteurs sélectionnés, (5) validation externe des modèles.

Modélisation QSPR de la constante de liaison des complexes liés à l'halogène.

L'ensemble de données comprenait 598 composés organiques pour lesquels le logarithme de la constante de stabilité du complexe 1:1 avec I 2 (logK BI2 ) a été mesuré dans l'hexane à 298K. Différents types de modèles de consensus correspondant à 4 scénarios possibles de génération de descripteurs à base de MA ont été comparés. Les meilleurs descripteurs, de type MA3, conduisent à des performances prédictives raisonnables à la fois dans la validation croisée et sur l'ensemble externe de 11 composés polyfonctionnels portant 2 ou 3 sites de liaison putatifs (Tableau 2 et Figure 2).

Modélisation QSPR de l'énergie libre de liaison hydrogène.

L'ensemble de données comprenait 3373 paires de complexes 1:1 formant une seule liaison hydrogène, pour lesquelles les mesures expérimentales de ΔG (kJ / mol) ont été reportées en conditions standard (dans du CCl 4 et à 298K). Les modèles consensus SVR et MLR, basés sur les descripteurs MA3 les plus performants, ont été utilisés pour la prédiction de l'ensemble de test externe de 629 complexes (R 2 = 0.65-0.74, RMSE = 3.2-3.8 (Figure 3)), mesurés dans différents solvants puis ramenés au CCl 4 , et l'ensemble d'essai de 12 complexes polyfonctionnels (Figure 4).

Modélisation QSPR et visualisation de la constante d'équilibre (logK T ) de différentes classes de transformations tautomères.

L'ensemble de données consistait en 695 équilibres tautomères attribués à 10 classes de transformation distinctes. Les deux modèles de classification prédisant le type de tautomérisation et les modèles de régression prédisant le logK T ont été obtenus avec les méthodes SVM et GTM. Les modèles impliquant des descripteurs locaux MA2 fonctionnent bien sur deux ensembles de tests externes: l'un contenant des transformations tautomères connues étudiées dans de nouvelles conditions (test set 1, Figure 5) et l'autre contenant de nouvelles transformations, au sens de leurs structures chimiques (test set 2, Figure 5).

Cet ensemble de données a été visualisé sur une carte bidimensionnelle construite en utilisant l'approche GTM (Figure 6). Cette carte sépare avec succès différentes classes de transformations tautomères et les équilibres différant soit par structure, soit par solvant.

Modélisation QSPR des propriétés cinétiques des réactions de cycloaddition.

Le jeu de données incluait 1849 réactions de types (4+2), ( 3+2 Les questions méthodologiques suivantes ont été considérées dans notre travail:

(1) une stratégie axée sur le processus de génération de descripteurs locaux; ( 2

Introduction

The requirements of modern knowledge-intensive fields of chemistry, such as drug development and chemical engineering, are constantly expanding, requiring more elaborated techniques and methods. Started by focusing on management of structural information of single molecular entities, present-day Chemoinformatics is developing the methods that follow the practical interest in predicting properties of compounds in condensed phase and in interaction with various partners, and not those of single molecules in vacuum. Consequently, tools and approaches maintaining single molecule study is not satisfying for multi-depended properties, i.e. the ones that depend on several parameters, or assemblage of molecules, bounded by a variety of intermolecular interactions. Exemplary objects are the host-guest complexation and molecular recognition, driven by weak interactions and, certainly, chemical reactions and chemical equilibria, thermodynamics and kinetics of which depends simultaneously on chemical structure and on reaction conditions. The mentioned processes represent the interactions incorporating several molecules, the structure of each of which has to be taken into thorough consideration. That, however, is not sufficient: for the case of intermolecular interactions (for instance), the interaction of two molecules, possessing several putative active sites, could result in different intermolecular complexes, that are not possible to differentiate if only the generic structures of the reagents are taken into account. The dynamics of the chemical process is thus characterized by the structural localities, directly involved into a chemical interaction and forming the active sites. Explicitly accentuated, the active sites define the dynamics of a particular chemical interaction. A generic task of prediction of the possibility of an interaction hence grow up into a new challenge of prediction of interactions with multiple putative active sites, for which one needs to predict which centers will interact and their interaction strength. This level requires a rigorous description of structural aspects of all the reagents, altogether with the accounting and explicit designation of the local regions signifying the active sites.

Another degree of freedom comes from the necessity of the reaction condition consideration. Indeed, small changes in solvent nature could drastically influence the property, in some cases being a determinative factor in the feasibility of the process. A simple solvent effect consideration may include a categorical assignment of a solvent to a certain group, such as polar/nonpolar or protic/aprotic. However, with regard to chemical reactions or intermolecular interactions, this simple approach could not already be sufficient. Indeed, the equilibrium or the rate constants remarkably depend on various factors included into the reaction conditions, such as solvent's polarity, solvent's H-acidity/basicity, temperature, pressure, etc. Thus, the new challenge addressed in this work is the necessity of the consideration and the description of the experimental conditions, including variety of solvent effects the property could depends on.

Chemical reaction data requires to be analyzed in terms of structural/conditional contents, its relationship with each other and with the property and data distribution in a chemical space. As an example, for the case of chemical reactions, the arrangement of the structural and reaction condition parts according with its influence on the property value, helps to understand the nature and the mechanism of the process. The need of such kind of analysis of complex data is increasing with the growth of the number of data constituents. A tool aimed at data analysis and visualization and used in this study is the Generative Topographic Mapping [START_REF] Bishop | GTM: The generative topographic mapping[END_REF][START_REF] Gaspar | GTM-Based QSAR Models and Their Applicability Domains[END_REF]8 (GTM). The method provides a visual, 2D-map projected data clustering and property distribution which is of great support for large data analysis. In this study, for the first time, it is used for complex chemical data modeling and analysis. Apart from common structural patterns identification, GTM helped to demonstrate the reaction condition influence, to estimate the quality of the data and the preferable way of its description. Moreover, the understanding of the underlying clustering principles contributes to revealing of the interrelation of data constituents and are helpful in determination of more suitable modeling approach and way of description.

One of the most demanded practical task, related to chemical processes, is the prediction of thermodynamic, kinetic or other parameters of a certain transformation. That could be achieved by the Quantitative Structure-Property Relationship (QSPR) modeling, one of the main tools of chemoinformatics, the goal of which is to provide an equation that relates an object (e.g. chemical reaction) with the value of the property of interest (e.g. rate constant). Various algorithms, so-called machine learning methods, have been developed for the QSPR modeling, each of which derive the equation in its own way. The mentioned GTM, first designed as the tool for data visualization, has been further extended for QSPR modeling, allowed a single GTM model to be used for the prediction of different, not necessarily related properties. The supervised methods such as Support Vector Machine (SVR) or Multiple Linear Regression (MLR), notwithstanding their inability in chemical data visualization, though could provide more accurate results of the prediction, since they are built specifically for a given property.

Thus, the challenge of this work is to expand QSPR methodology for problems, where the working hypothesis of a "constant" environment does no longer apply: interactions of both covalent and non-covalent nature, with different partners in multiple solvents, at various temperatures. These processes are characterized by local interactions, that imply the representation of both, structural and conditional aspects, along with the explicit designation of the dynamics of the process. The complexity of the tasks is raising form intermolecular interactions to chemical reactions.

The work is divided into two parts, where the first one describes the general methodology and the second is devoted to its practical application for different chemical objects. The first chapter gives an overview of the existing local descriptors, providing the structural characterization of a process, and the descriptors encompassing the solvent effects. The purpose of that section is to select an appropriate way of chemical data representation. As it will be discussed in the next section, one of the most convenient types of description for complex chemical processes is the ISIDA Marked Atom-based (MA) and the ISIDA Condensed Graph of Reaction-based (CGR) fragment descriptors, representing substructures of a molecular structure. The second chapter describes the general practices of QSPR modeling and the details and particularities of the machine learning methods, used during the study: SVM, MLR and GTM. The second part of the thesis, devoted to the practical application, consists of five projects. The part of intermolecular interaction modeling starts with halogen bonding, where the data represent a set of single molecules measured in unified conditions. The topic continues with hydrogen bonding interaction, for which the model predicting the strength of intermolecularly-bonded complexes of different donors and acceptors was built. The work continues with even more challenging modeling of chemical reactions. The section starts with the project of tautomeric equilibria modeling and visualization, accounting for the impact of reaction condition changes. The section continues with an exhaustive modeling and visualization of kinetic parameters of reactions of cycloaddition. The final chapter is dedicated to the modeling of a large set of S N 1 reactions, where both approaches of structure representation (CGR-based and MA-based) are employed. Table 1 illustrates the overview of the application part providing the class of chemical processes, property and the descriptors and tools that have been employed during the study. i. Glavatskikh, M., Madzhidov, T., Solov'ev, V., Marcou, G., Horvath, D., Graton, J., ... & Varnek, A. (2016). Predictive Models for Halogen-bond Basicity of Binding Sites of Polyfunctional Molecules. Molecular informatics, 35(2), 70-80. ii.

Glavatskikh, M., Madzhidov, T., Solov'ev, V., Marcou, G., Horvath, D., & Varnek, A. (2016). Predictive models for the free energy of hydrogen bonded complexes with single and cooperative hydrogen bonds. Molecular informatics, 35 (11-12), 629-638.

Chapter 2

Molecular descriptors for interacting chemical entities

The Quantitative Structure-Property Relationship (QSPR) presupposes a molecular structure to be encoded in a way that is, first, adapted for the machine learning algorithm and, second, convenient for the evaluation of the corresponding structure-property relationship. The attributes used for the description of a molecule, called descriptors, should be chosen in accordance with the task, taking into account the nature of the process, its driving force and the factors that could affect the process. Regarding to the modeled processes, a simplified categorization could include global processes, referred to the whole structure, for which, among the scope of possible local interactions, the property-defined ones could not be specified (solubility). These properties are thus could be considered as depending on the structure as a whole. Local processes are defined primarily by local interactions of the active centers (intermolecular binding). Thus, for local processes the interaction centers could be pointed out.

This chapter gives an outlook of the variety of descriptors applicable for the QSPR modeling of local processes (local descriptors). The description of a chemical process includes the structural and the experimental condition parts. Correspondingly, local descriptors are referred to chemical structure, and discussed in section 2.1, whereas the parameters, by which the reaction conditions could be taken into account, are reviewed in section 2.2. A comprehensive overview of this field is not a goal of the present chapter and the attention will be paid to the most common, widely-used descriptors conforming the definition of 'local descriptors' i.e. bearing an information about certain atoms or group of atoms. That includes: substituent constants, quantum-chemical descriptors, electrotopological indices and ISIDA fragment descriptors. A full comprehensive review of major types of global and local descriptors, used in chemoinformatics, are given in works and books of R.Todeschini and V.Consonni 9-10 .

Local descriptors for chemical structure representation

Substituent constants

Substituent constants could be proclaimed as the first attempt in classification of local effects of certain structural groups. A pioneer work belongs to Hammet 11 , treating the electronic effect of substituents on the rate and equilibrium constants of organic reactions, and Taft [START_REF] Taft | The evaluation of inductive and resonance effects in reactivity. II. Thermodynamic properties of hydrogenation of non-conjugated olefins, aldehydes and ketones[END_REF] , applying similar approach for derivation of series of constants, differentiated by nature of the contributed electronic effect.

Hammet constants

An American physical chemist Louis Hammett noted that a particular substituent on the aromatic ring of benzoic acid would affect its acidity in a similar manner as it would affect other aromatic structures. For instance, a para-nitro group would affect the value of the dissociation of benzoic acid in a manner similar to that of salicylic acid. That was the beginning of the concept of substituent constants. The well-known Hammet constants are derived from the dissociation constants ratio of benzoic acid (K 0 ) and a corresponding substituted benzoic acid:

σ = log K K 0 ⁄ (1)
Because of the drastic dependence of the dissociation constants upon temperature and the nature of solvent, the σ-constants are specifically given for water solution and at the temperature of 25° C.

The magnitude of the electronic effect caused by the substituent is influenced by its position to the carboxylic group. In this way, the σ-constant for para-position will mainly describe the electronic influence by means of resonance effect, while the σ-constant for ortho-position will describe both fluctuations via σ and π -bonds. Since the strength of the effects varies depending on the position of the substituent in the ring, the meta, ortho and para constants, σ m, σo and σp, are distinguished. . If the ratio K K 0 ⁄ is more than one, i.e., the substituent leads to an increased acidity of the benzoic acid, σ is positive and the substituent is considered to be an electron-withdrawing group, if the ratio is less than one, the substituent is electrondonating and σ will be negative. Hammett substituent constants are referred to hydrogen and σ H is thus equal to zero.

Despite of the empirical derivation from the benzoic acids equilibrium, substituent constants can be successfully applied for the prediction of the variety of families of reactions in solution, such as electrophilicity of substituted benzoic esters, the nucleophilicity of anilines, and the solvolysis of benzyl halides [START_REF]A guidebook to mechanism in organic chemistry[END_REF] . Hammet constants are important constituents in the field of QSPR modeling and were applied for the modeling of protein-ligand interactions [START_REF] Sung | Functional analogue reaction systems of the DMSO reductase isoenzyme family: Probable mechanism of S-oxide reduction in oxo transfer reactions mediated by bis(dithiolene)-tungsten(IV,VI) complexes[END_REF] , interactions with enzymes [START_REF] Hansch | Molecular graphics and QSAR in the study of enzyme-ligand interactions. On the definition of bioreceptors[END_REF] , antitumor and antimalarial activity [START_REF] Selassie | Comparison of triazines as inhibitors of L1210 dihydrofolate reductase and of L1210 cells sensitive and resistant to methotrexate[END_REF][START_REF] Wipf | Synthesis and biological evaluation of the first pentafluorosulfanyl analogs of mefloquine[END_REF] as well as toxicology and mutagenicity [START_REF] Hansch | A survey of Hammett substituent constants and resonance and field parameters[END_REF] .

Inductive constants

The electronic constants devised by Hammett reflects three types of electronic influences:

• resonance (mesomeric) effect

• inductive effect: electrostatic influence of a group which is transmitted primarily by polarization through a chain of neighboring atoms

• field effect: electrical influence of a substituent transmitted through space

The last two are hard to distinguish and usually they are considered to be a unified composite inductive effect and are treated together. Thus, because of the complexity and unified nature of the overall electronic constants, the establishing of the way by which the substituent influences on the reaction rate and equilibrium is very important, because some chemical reactions are driven either by the resonance or by the inductive effect. The approaches of quantitative evaluation of a pure inductive effect were first devised by Taft and Ingold [START_REF] Taft | The general applicability of a fixed scale of inductive effects. II. Inductive effects of dipolar substituents in the reactivities of m-and p-substituted derivatives of benzene 1,2[END_REF] and then proceeded in works of Roberts and Moreland 20 , Holtz-Stock 21 , Siegel-Kormany [START_REF] Siegel | Quantitative relationships in the reactions of trans-4-X-cyclohexanecarboxylic acids and their methyl esters[END_REF] and others.

Taft inductive constant

In 1930 Ingold proposed the idea of measurement of an inductive effect through a ratio of dissociation constants rates of acid and base hydrolysis of acetic acid esters. Developing the Ingold's idea, Taft [START_REF] Taft | The general applicability of a fixed scale of inductive effects. II. Inductive effects of dipolar substituents in the reactivities of m-and p-substituted derivatives of benzene 1,2[END_REF] derived series of inductive substituent constants (σ * constants) estimating quantitatively an inductive effect and defined as :

𝜎 * = 1 2.48 ⁄ × [log( 𝑘 𝑥 𝑘 𝑀𝑒 ) 𝑏 -log( 𝑘 𝑥 𝑘 𝑀𝑒 ) 𝑎 ] (2)
where the indexes A and B refer to acid and base hydrolysis. The factor 2.48 is introduced to make the 𝜎 * values comparable in magnitude to the widely used Hammett constants. A positive 𝜎 * value indicates that the group is electron-withdrawing relative to methyl, while a negative value indicates electron contribution. In acid and base series of the reactions, the steric and resonance effects can be considered to be the same: the transition state of both mechanisms passes through tetrahedral intermediate (Figure 1). The acidic intermediate will differ from the one of base-catalyzed by one proton, which can not affect the steric factor significantly. In case of the resonance influence possibility, it will also be involved in both intermediates and the effect should be nearly the same for both mechanism.

Roberts-Moreland inductive constant

The Robert-Moreland constants [START_REF] Roberts | σ-constants of the carbethoxyl and hydroxyl groups[END_REF] derived from the measurement of the dissociation constants for a series of 4-substituted bicyclo-[2.2.2]-octane-l-carboxylic acids (Fig. 2). This molecule has no unsaturation, hence, the transmission of electrical effects of substituents through the ring by resonance is not possible and the substituent can induce the inductive effect only. Moreover, the chosen reference compound is free from conformational effects and no steric effect is observed, as the substituent and the active site are not in close proximity to each other. The dissociation constant is measured in 50% ethanol at 25°C.

The Roberts-Moreland inductive constant measured in 50% ethanol at 25°C, and defined as:

𝜎 = 1 1.464 ( 𝑙𝑜𝑔𝐾 -𝑙𝑜𝑔𝐾 0 ) (3)
Where 𝐾 0 is the dissociation constant of unsubstituted bicyclo-[2.2.2]-octane-l-carboxylic acid. The coefficient of 1.464 is given to refer the scale to the Hammett equation.

Resonance (mesomeric) constants

Taft resonance constant

The application of the well-known Hammett sigma constants referred to meta-and parasubstitutions sometimes can be limited by the fact that some reactions are mostly driven by resonance effect, which is implicitly included into the para-substitution constant and not for the meta-one. Thus, the resonance contribution for the series of para-substituted benzene derivatives can be simply expressed through the subtraction of the pure inductive contribution from the Hammett sigma constant. That was done by Taft 12 who proposed the first scale of a resonance effect of the given series of compounds:

𝜎 𝑟 = 𝜎 -𝜎 𝑖 (4)
Where 𝜎 𝑖 is the Taft inductive constant (see 2.1.1.2). The resonance constants express the influence of the π-bonded electrons of the substituent to the benzene ring due to resonance fluctuation. As a measure of withdrawing of the electronic charge, the values of the 𝜎 𝑟 are negative for ortho-and para-groups and positive for the meta-groups. Later on, Taft 23 also proposed the estimation of 𝜎 𝑟 by means of [START_REF] Taft | The general applicability of a fixed scale of inductive effects. II. Inductive effects of dipolar substituents in the reactivities of m-and p-substituted derivatives of benzene 1,2[END_REF] F-NMR spectroscopy, where the 19 F-chemical shift indicates the resonance interaction between the para-substituent and fluorobenzene system.

It should be noticed, that the 𝜎 𝑟 resonance scale is only suitable for benzene derivatives as the resonance effect in general has a great variability upon the reaction type. Although, the scale can provide a general qualitative estimation of resonance ability of a certain substituent.

Swain-Lupton approach

The other approach of separate estimation of the inductive and resonance effects was proposed by Swain and Lupton. The idea of the authors is that the Hammett constant, which included inductive along with the resonance effect, can be represented as a linear combination of both contributions with the corresponding parameters:

𝜎 = 𝑓𝐹 + 𝑟𝑅 (5)
where the polar component (𝐹) is calculated from the 𝜎 𝑚 and 𝜎 𝑝 Hammett constants: 𝐹 = 𝑏 0 + 𝑏 1 𝜎 𝑚 + 𝑏 2 𝜎 𝑝 (the coefficients are evaluated by least square regression using pK a values of bicyclo-[2.2.2]-octane-l-carboxylic acid), and the resonance component 𝑅 is estimated as 𝜎 𝑝 -0.921𝐹.

The main assumption hence that the substituent in para-position induce the main resonance perturbation. The corresponding 𝐹 and 𝑅 parameters were initially calculated for 43 substituents and then further expanded to a few hundreds.

Steric constants

Taft steric constant

The first steric constant 𝐸𝑠 was defined empirically by Taft as the extension of Hammett equation [START_REF] Hammett | The effect of structure upon the reactions of organic compounds. Benzene derivatives[END_REF] . 𝐸𝑠 is a measurement of the steric effect caused by the group X and influenced the acid-catalyzed hydrolytic rate of esters of substituted acetic acids: 𝐸𝑠 = log (𝑘 𝑋 ) 𝐴 -log (𝑘 𝐻 ) 𝐴 (6) where 𝑘 𝑋 and 𝑘 𝐻 are the rates of substituted and unsubstituted acetic acids esters hydrolysis. This scale is based on the assumption that the corresponding rates are influenced mainly by steric effects and no polar interruption is introduced. The bulkier the substituent, the more negative the 𝐸𝑠 constant value is.

The 𝐸𝑠 scale succeeded in reproducing of steric effect, giving, at least, qualitative approximation for the measured substituent effect. Later on, more unified and revised scales have been proposed: Hancock 24 corrected the 𝐸𝑠 parameter with the inclusion of the hyperconjugation influence of α-hydrogens, Palm [START_REF] Palm | Fundamentals of Quantitative Theory of Organic Reactions[END_REF] enlarged the latter with the C-C and C-H hyperconjugation effect corrections, Dubois [START_REF] Hirota | A new steric substituent constant Ωs based on molecular mechanics calculations[END_REF] proposed a modified scale defined on the basis of more unified reactions and over wider range of substituents.

Charton steric constant

Charton found that Taft's steric constant is linearly dependent on the van der Waals radius of substituent, which led to the developments of Charton's steric parameter. The constant is defined as the difference of the corresponding van der Waals radius of substituent X and hydrogen atom radius:

𝑣 𝑥 = 𝑅 𝑣𝑑𝑤(𝑋) -𝑅 𝑣𝑑𝑤(𝐻) (7) υ is defined as the difference between the van der Waals radii of H and substituent X.

Charton's steric parameter related to the van der Waals radius of any symmetrical substituent (H, Cl, CN) or to the minimum width of asymmetrical ones (CH3, CMe3). Charton also defined the minimum and maximum van der Waals radius in order to take into account the possibility of conformation of a group thus seeking for the repulsive effect minimization, the average of which well correlated to Taft steric constant.

Quantum-chemical descriptors

Atomic charges

According to the classical chemical theory, the driving force of all chemical reactions is either of the electrostatic or of the orbital-control driven nature. Thus, charges are responsible for the whole variety of electrostatic-driven processes. It has been shown that local electron densities or charges are essential parameters in description and interpretation of the mechanism of chemical reactions and physico-chemical properties [START_REF] Karelson | Quantum-chemical descriptors in QSAR/QSPR studies[END_REF] . That is the reason of wide usage of charge-based descriptors in QSPR modeling of different physico-chemical properties, chemical reactions and weak intra-or intermolecular interactions.

Most common schemes for atomic charges derivation are based on the population analysis of the wave function obtained by quantum-chemical calculation. Several schemes for the analysis of the wave function have been proposed. The most common and utilized are Mulliken [START_REF] Mulliken | Electronic Population Analysis on LCAO-MO Molecular Wave Functions[END_REF] and Löwdin [START_REF] Löwdin | On the Nonorthogonality Problem*[END_REF] atomic charges, those based on natural bond orbital theory [START_REF] Foster | Natural hybrid orbitals[END_REF] (NBO), the Bader AIM theory [START_REF] Bader | A quantum theory of molecular structure and its applications[END_REF] and the ones fitting the point charges such as to produce an intrinsic electrostatic potential calculated from the wave function [START_REF] Breneman | Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis[END_REF] . The diversity of the calculation methods is a consequence of the fact that none of the values obtained by any of the methods corresponds to a directly experimentally measurable quantity. That should be mentioned, however, that partial charges could be obtained by empirical methods, such as Gasteiger-Marsilli 33 .

Atomic charges have been used as static chemical reactivity indices. One of the most commonly used nondirectional indices are net atomic charges, which can be obtained by subtracting the number of valence electron belonging to the atom from the total electron density on the atom. As a global version of charge-based descriptors, the most positive and the most negative net atomic charges and the average absolute atomic charge are often used 34- 35 .

Atomic charges have been successfully used as local descriptors for QSPR modeling of different physico-chemical properties such as partition octanol-air coefficient [START_REF] Chen | Quantitative structure-property relationships for octanol-air partition coefficients of PCDD/Fs[END_REF] , adsorption coefficient [START_REF] Baker | Evaluation of estimation methods for organic carbon normalized sorption coefficients[END_REF] , dopamine and benzodiazepine agonists [START_REF] Bauknecht | Locating biologically active compounds in medium-sized heterogeneous datasets by topological autocorrelation vectors: Dopamine and benzodiazepine agonists[END_REF] , acid dissociation constant (pKa) [START_REF] Dixon | Estimation of pKa for organic oxyacids using calculated atomic charges[END_REF][START_REF] Svobodova | Predicting pK(a) values of substituted phenols from atomic charges: comparison of different quantum mechanical methods and charge distribution schemes[END_REF] and hydrogen-bong strength prediction [START_REF] Gancia | Theoretical hydrogen bonding parameters for drug design[END_REF] .

Electrophilic and nucleophilic frontier electron densities

One of the most powerful tool for chemical reactivity interpretation is the frontier molecular orbitals theory (FMO), developed by Kenichi Fukui in 1950's. The theory is based on the consideration of the frontier molecular orbitals, correspondingly, the highest occupied and the lowest unoccupied molecular orbitals (HOMO and LUMO), as the ones mainly responsible for molecule's reactivity. Thus, the frontier orbital theory predicts the site of the lowest unoccupied orbital localization to be an electrophilic region, similarly, the site where the highest occupied orbital is localized is a nucleophilic region.

The theory gave rise to many different global and local descriptors which are widely usable due to its high information content, wide applicability and easiness of calculation. The most common local FMO descriptors are based on the atomic contribution to HOMO or LUMO.

Thus, an electrophilic frontier electron density 𝐹 𝑎 𝐸 indicates how easy the atom a interacts with an electrophile. Opposite, a nucleophilic frontier electron density 𝐹 𝑎 𝑁 is a measure of the atom a to be exposed for the nucleophilic attack. These descriptors are defined as:

𝐹 𝑎 𝐸 = ∑(𝐶 ℎ𝑜𝑚𝑜 ,𝑗) 2 |𝐸 ℎ𝑜𝑚𝑜 | and 𝐹 𝑎 𝑁 = ∑(𝐶 𝑙𝑢𝑚𝑜 ,𝑗) 2 |𝐸 𝑙𝑢𝑚𝑜 | ( 8
)
where 𝐶 ℎ𝑜𝑚𝑜 , 𝑗 and 𝐶 𝑙𝑢𝑚𝑜 , 𝑗 are the coefficients of contributions of the j-atomic orbital of the atom a to HOMO and LUMO, and 𝐸 ℎ𝑜𝑚𝑜 𝐸 𝑙𝑢𝑚𝑜 are the energies of the corresponding orbitals. The FMO descriptors perform better when HOMO and LUMO are well separated in energy and the reaction is fully controlled by the frontier orbitals (which, for example, is not the case of aromatic ring system). The examples of the application of the electrophilic and nucleophilic frontier electron densities are: modeling of mutagenicity [START_REF] Tuppurainen | About the mutagenicity of chlorine-substituted furanones and halopropenals -a QSAR study using molecular-orbital indexes[END_REF] , antioxidant activity [START_REF] Sarkar | A QSAR study of radical scavenging antioxidant activity of a series of flavonoids using DFT based quantum chemical descriptors -the importance of group frontier electron density[END_REF] , adsorption of organic compounds on soils [START_REF] Liu | QSAR analysis of soil sorption coefficients for polar organic chemicals: Substituted anilines and phenols[END_REF] and porphines and chlorins reactivities [START_REF] Knop | The reactivities of porphin, chlorin, bacteriochlorin, and phlorin. Electron densities, free valences and frontier orbital densities[END_REF] .

Electrophilic, nucleophilic and radical superdelocalizabilities

Along with the electrophilic and nucleophilic frontier electron densities, another type of descriptors derived from the Fukui's theory are superdelocalizability indices, which can be defined as the contribution of the atom a to the stabilization energy in the formation of a charge-transfer complex or the ability to form bonds through charge transfer. Thus, the electrophilic superdelocalizability (𝑆 𝑎 𝐸 ) describes the interaction with the electrophilic center and the nucleophilic superdelocalizability (𝑆 𝑎 𝑁 ) describes the interaction with the nucleophilic center:

𝑆 𝑎 𝐸 = 2 × ∑ ∑ 𝑐 𝑖,𝑗 2 |𝐸 𝑖 | 𝑁𝑜𝑐𝑐 𝑖 𝑎𝑛𝑑 𝑆 𝑎 𝑁 = 2 × ∑ ∑ 𝑐 𝑖,𝑗 2 |𝐸 𝑖 | 𝑁𝑚𝑜 𝑁𝑜𝑐𝑐+1 ( 9
)
where 𝑐 𝑖,𝑗 are the coefficients of the contribution of the j's atomic orbital to the i's molecular orbital summed over all occupied/unoccupied molecular orbitals, 𝐸 𝑖 is the energy of the i's molecular orbital.

These are useful parameters to characterize molecular interactions and to compare corresponding atoms in different molecules. Electrophilic and nucleophilic superdelocalizabilities themselves are local descriptors that describe atom in a molecule, in parallel, there are important global descriptors, based on the atomic superdelocalizabilities such as maximum, total and average superdelocalizabilities. Superdelocalizabilities are socalled dynamic reactivity indices, referring to the transition states of the reactions [START_REF] Karelson | Quantum-Chemical Descriptors in QSAR/QSPR Studies[END_REF] , while the static indices (e.g. charges) describe isolated molecules in their ground state.

Superdelocalizability indices have been used as descriptors in different works devoted to the modeling of physico-chemical parameters or reactivities. Some of the examples are modeling of the acute toxicity of the substituted benzenes [START_REF] Karabunarliev | Quantum-chemical descriptors for estimating the acute toxicity of substituted benzenes to the guppy (Poecilia reticulata) and fathead minnow (Pimephales promelas)[END_REF] , modeling of series of benzodioxanes as alpha-1-adrenergic antagonists [START_REF] Venturelli | Molecular modeling and quantitative structure activity relationship analysis using theoretical descriptors of 1,4-benzodioxan (wb-4101) related-compounds alpha-1-adrenergic antagonists[END_REF] , of carbonic anhydrase inhibitors [START_REF] Supuran | A quantitative structure-activity relationship study of positively charged sulfonamide inhibitors[END_REF] and of the permeability coefficient of aminobenzoates [START_REF] Fu | Study on percutaneous rate of P-aminobenzoates using molecular orbital method[END_REF] .

Radical superdelocalizability is defined as the sum of the electrophilic and nucleophilic superdelocalizabilities:

𝑆 𝑎 𝑅 = 2 × ∑ ∑ 𝑐 𝑖,𝑗 2 |𝐸 𝑖 | 𝑁𝑜𝑐𝑐 𝑖 + 2 × ∑ ∑ 𝑐 𝑖,𝑗 2 |𝐸 𝑖 | 𝑁𝑚𝑜 𝑁𝑜𝑐𝑐+1 ( 10
)
Radical superdelocalizability refers to radical attack and have been used as an important descriptor in the modeling of toxicity of halogenated aliphatic compounds [START_REF] Chroust | Quantitative structure-activity relationships for toxicity and genotoxicity of halogenated aliphatic compounds: Wing spot test of Drosophila melanogaster[END_REF] , reactions of hydroxyl radicals with nucleic acids [START_REF] Masuda | Reactions of hydroxyl radicals with nucleic acid bases and the related compounds in gamma-irradiated aqueous solution[END_REF] , reactions of hydroxylations of aromatic compounds [START_REF] Omori | Relation between electronic structure and hydroxylation of aromatic compounds by microorganisms[END_REF] and carcinogenicity of polycyclic hydrocarbons [START_REF] Pullman | Electronic aspects of the interactions between the carcinogens and possible cellular sites of their activity[END_REF] .

Atomic polarizability

Among common local quantum-chemical descriptors, the polarizability indices occupy an important place. In general, atomic polarizability is the polarization effect at atomic level, where dipole moment µ 𝑖𝑛𝑑,𝑖 is induced on the ith atom:

µ 𝑖𝑛𝑑,𝑖 = 𝑎 𝑖 × 𝐸 𝑖 (11) where 𝐸 𝑖 is the electric field at the ith atom and 𝑎 𝑖 is the corresponding atomic polarizability tensor.

Several methods have been proposed for the atomic polarizability calculation. One of the first was developed by Kang 55 in which the atomic polarizabilities were obtained from the experimental polarizabilities of homologous molecules. The method of Miller [START_REF] Miller | Additivity methods in molecular polarizability[END_REF] allows to calculate so-called atomic hybrid polarizabilities which take into account the hybridization of the atom. These atomic hybrid polarizabilities can be combined to generate bond polarizabilities and the average molecular polarizability. The method developed by [START_REF] No | An empirical-method to calculate average molecular polarizabilities from the dependence of effective atomic polarizabilities on net atomic charge[END_REF] proposes to calculate the effective atomic polarizabilities as functions of net atomic charges.

In addition to the simple atomic polarizabilities, the common descriptors of this family are atom-atom polarizability and self-atomic polarizabilities. Atom-atom polarizability is an index of chemical reactivity, denoted as 𝜋 𝑎𝑏 and calculated from the perturbation theory as: where i and j run over the molecular orbitals and µ and v run over the atomic orbitals, 𝑐 𝑖µ,𝑎 denotes the i-th molecular orbital coefficient for atomic orbital μ located on atom a.

𝜋 𝑎𝑏 = 4 ∑ ∑ ∑ ∑ 𝑐 𝑖µ
The self-atom polarizability is analogously defined as (13) Polarizability indices have been successfully applied for the calculation of the conjugation energies [START_REF] Brown | Conjugation energy[END_REF] , nuclear spin-spin coupling constants [START_REF] Joela | Theoretical nuclear spin-spin coupling constants using atomatom polarizabilities. 13C-H and H-H′ coupling constants of some [2.2.1bicyclic compounds using the INDO approximation[END_REF] , treatment of induction effects in molecular mechanics simulations [START_REF] Soteras | Derivation of distributed models of atomic polarizability for molecular Simulations[END_REF] and carcinogenicity of nitroso-compounds [START_REF] Helguera | Quantitative structure carcinogenicity relationship for detecting structural alerts in nitroso-compounds. Species: Rat; Sex: Male; Route of administration: Water[END_REF] .

𝜋 𝑎𝑏 = 4 ∑ ∑ ∑ ∑ 𝑐 2 𝑖µ,𝑎 𝑐 2 𝑗𝑣,𝑏 𝐸 𝑖 -𝐸 𝑗 𝑣 µ 𝑁𝑜𝑐𝑐 𝑗 𝑁𝑜𝑐𝑐 𝑖

2.1.2.5

TAE descriptors based on Bader's quantum theory of atoms in molecules

The theory of Atoms in Molecules (AIM) was developed by Bader [START_REF] Bader | A topological theory of molecular-structure[END_REF] and remains to be commonly used and applicable methods for the calculation of atomic and different molecular properties and study of molecular interactions. The theory is based on the properties of the observable charge distribution of a molecular system and provides a unique mapping between the topological elements of a molecular charge distribution and the structural elements, atoms and bonds, underlying the notion of molecular structure [START_REF] Matta | An Introduction to the Quantum Theory of Atoms in Molecules[END_REF] . Central to this theory is the identification of an atom with a particular region of real space as determined by a fundamental topological property of a charge distribution. By appealing to quantum mechanics one finds that the atoms so defined possess a unique set of properties and behave as closed physical system. In particular, the theory shows that the average value of every mechanical property (a property whose associated operator can be expressed in terms of the coordinate and/or momentum operators) of some system can be expressed as a sum of corresponding atomic contributions. The total energy of a crystal, for example, is equal to the sum of the energies of the atoms in the crystal where each atom is a well-defined object in real space. An important point of the theory that properties attributed to atoms and functional groups are transferable from one molecule to another [START_REF] Bader | Transferability of atomic properties and the theorem of Hohenberg and Kohn[END_REF] . The most applicable and practically used characteristics coming from AIM theory are bond critical point properties and atomic properties.

Bond critical points (BCP) are saddle points in electron density distribution in the region between bonded atoms having two negative and one positive eigenvalue of hessian. Several BCP properties have been shown to be correlated with experimental molecular properties [START_REF] Becke | The quantum theory of atoms in molecules: from solid state to DNA and drug design[END_REF] . For example, the electron density at the BCP correlates with the bond energies, and hence provides a measure of bond order 66 , the potential energy density at the BCP has been shown to be highly correlated with hydrogen bond energies [START_REF] Espinosa | Hydrogen bond strengths revealed by topological analyses of experimentally observed electron densities[END_REF] and theoretically computed proton shielding [START_REF] Buttingsrud | An investigation of descriptors based on the critical points in the electron density by building quantitative structure-property relationships for proton chemical shifts[END_REF] .

Atomic properties have been used to recover and directly predict several additive atomic and group contributions to molecular properties, including, for example, heats of formation 69 magnetic susceptibility [START_REF] Keith | Calculation of magnetic response properties using atoms in molecules[END_REF] , molecular volumes [START_REF] Bader | Properties of atoms in molecules -atomic volumes[END_REF] , dipole moment 72 , polarizability [START_REF] Gough | Analysis of molecular polarizabilities and polarizability derivatives in H-2, N-2, F-2, CO, and HF, with the theory of atoms in molecules[END_REF][START_REF] Bader | Properties of atoms in molecules -additivity and transferability of group polarizabilities[END_REF] and many others. Atomic properties have also been used build QSPR models predicting several experimental properties including, for example, the pKa of carboxylic acids, anilines and phenols [START_REF] Chaudry | Estimation of pK(a) using quantum topological molecular similarity descriptors: Application to carboxylic acids, anilines and phenols[END_REF] a wide array of biological and physicochemical properties of the amino acids, and the effects of mutation on protein stability [START_REF] Matta | Atoms-in-molecules study of the genetically encoded amino acids. III. Bond and atomic properties and their correlations with experiment including mutation-induced changes in protein stability and genetic coding[END_REF] , NMR spin-spin coupling constants of aromatic compounds from the electron delocalization indices [START_REF] Castillo | Fluorine-fluorine spin-spin coupling constants in aromatic compounds: Correlations with the delocalization index and with the internuclear separation[END_REF] .

However, the properties of atoms and bonds derived from the QTAIM and based on quantummechanical calculation require significant computational costs. In order to overcome the problem, Breneman [START_REF] Breneman | Electrondensity modeling of large systems using the transferable atom equivalent method[END_REF] introduced the concept of 'transferable atom equivalents' (TAE)atom-based electron density fragments obtained using the AIM approach. The underlying concepts for the TAE method is the additivity of atomic properties in Bader's theory and the transferability of topological atoms. A TAE is a mononuclear atomic region of space filled with electron density delimited by zero-flux surfaces in the gradient vector field of the electron density extracted from a parent molecule at its zero-flux surface. Extracted atoms representing a large number of differing combinations of elements, atom types, and immediate electronic environments, are stored in a computerized database. A program called RECON is then used to assemble the electron densities (and other properties) of a target large molecule by matching the appropriate zero-flux surfaces of different TAEs recalled from the database [START_REF] Matta | Electron-density descriptors as predictors in quantitative structure-activity/property relationships and drug design[END_REF] . Once the target molecule is reconstructed by the automated merging of TAEs, the molecular descriptors are then calculated by arithmetic or vector sums of the properties of the composing TAEs. The reconstructed descriptors include, for example, the total molecular energy, the total molecular volume, the electrostatic potential, the molecular dipole moment, and the Fukui functions. The TAE approach, being a useful tool in QSAR/QSPR modeling, has been proved to be relevant for modeling of protein-ligand binding affinity [START_REF] Ryan | Application of novel refined distance-dependent descriptors, Transferable Atom Equivalent (TAE) techniques and machine-learning methods in predicting protein-ligand binding affinity[END_REF] , Mu-opioid receptor affinity [START_REF] Shen | Modeling the Mu-opiod receptor affinity of synthetic 8-aminocyclazocine analogues using TAE, PEST and PAD descriptors and machine-learning methods[END_REF] , for high-throughput screening [START_REF] Sukumar | Virtual high-throughput screening of large datasets using TAE/RECON descriptors[END_REF] and molecular surface autocorrelation analysis [START_REF] Breneman | New developments in PEST shape/property hybrid descriptors[END_REF] .

Conceptual Density Functional Theory Indices

Among others commonly used local descriptors, Fukui Functions are one of the most popular in describing molecule's site selectivity and chemical reactivity. Fukui Functions find their origin within Conceptual Density Functional Theory (Conceptual DFT) and are defined as:

𝑓(𝑟) = 𝑑𝒑(𝑟) 𝑑𝑵(𝑟) = 𝑑µ 𝑑𝝂(𝑟) ( 14
)
where 𝒑(𝑟) is the electron density at a point r, 𝑵(𝑟) is a total number of electrons of the system at a given external potential 𝝂(𝑟). Besides, the Fukui function corresponds to the first derivative of the electronic chemical potential µ with respect to the external potential 𝝂(𝑟)

for a given number of electrons. Depending on the nature of the electron transfer, the Fukui function for removal of an electron from the molecule, called the Fukui function for electrophilic attack, is labeled as 𝑓 -, and the Fukui function for addition of an electron to the molecule, called the Fukui function for nucleophilic attack, is labeled as 𝑓 + are distinguished:

𝑓 -(𝑟) = 𝒑(𝑟) 𝑁 -𝒑(𝑟) 𝑁-1 (15)
𝑓 + (𝑟) = 𝒑(𝑟) 𝑁+1 -𝒑(𝑟) 𝑁 ) (16) where 𝒑(𝑟) 𝑁 is the electron density at a point r for the molecule possessing N electrons (N corresponds to neutral molecule). Thus, 𝑓 -is large in the regions of space where a given molecule readily donates electrons and 𝑓 + is large in the regions where a molecule accepts electrons. A reaction thus is likely to occur between regions (or atoms) where 𝑓 -is large in one molecule and 𝑓 + is large in another reacting molecule.

The evaluation of the Fukui function values is not straightforward and number of methods and algorithms have been developed in order to ease the calculation. Thus, Yang and Mortier [START_REF] Yang | The use of global and local molecularparameters for the analysis of the gas-phase basicity of amines[END_REF] proposed three different condensed forms of f (r), based on atomic charges of N, N + 1, and N -1 electron systems, Nalewajski [START_REF] Nalewajski | Qualitative charge stability and fukui function-analysis of a model system with implications for chemical-reactivity[END_REF] has studied the f (r) indices in respect of Bader's 'atom in molecules' (AIM) theory, Komorowski et al [START_REF] Komorowski | Electronegativity and hardness of chemical groups[END_REF] proposed the atomic and group resolution of f (r) indices based on semiempirical method. The most popular method, proposed by Yang and Mortier 84 , based on the condensation of the Fukui functions to atomic resolution:

𝑓 -(𝑟) = 𝒒 𝒂 (𝑁) -𝒒 𝒂 (𝑁 -1) (17)
𝑓 + (𝑟) = 𝒒 𝒂 (𝑁 + 1) -𝒒 𝒂 (𝑁) (18) where 𝒒 𝒂 is the charge on an atom a for a molecule having N electrons. The method has a simple procedure to calculate the atomic condensed Fukui function indices using a charge partitioning schemes, e.g. Natural, Mulliken or Hirschfeld Population Analysis. However, despite of the possibility of using any type of atomic charges, it was shown that Hirschfeld charges are likely the most accurate [START_REF] Hocquet | Intramolecular interactions along the reaction path of keto-enol tautomerism: Fukui functions as local softnesses and charges as local hardnesses[END_REF] for Fukui indices calculation. Thus, ranking of atoms within a molecule in terms of condensed Fukui functions enable the identification of preferential sites of reactions. Nevertheless, one should remember that the Fukui functions have a poor performance in handling the hard-hard interactions, but they are the good descriptors for the soft-soft interactions known to be frontier controlled.

Characteristic examples of QSRP modeling with Fukui functions included modeling of ketoenol tautomerism [START_REF] Hocquet | Intramolecular interactions along the reaction path of keto-enol tautomerism: Fukui functions as local softnesses and charges as local hardnesses[END_REF] , local reactivities during electrophilic, nucleophilic and radical attacks [START_REF] Kolandaivel | Study of atomic and condensed atomic indices for reactive sites of molecules[END_REF] and reactivity for protonation reactions [START_REF] Melin | Is the Fukui Function a Right Descriptor of Hard-Hard Interactions[END_REF] .

The topic of hard-soft interactions and the derived local chemical reactivity implies a logical continuation into local softness/hardness introduction. The concepts of a local softness and hardness, as goes from the name, is the exertion of the principles of softness and hardness in local sense so as to explain the response of a chemical system to different kinds of reagents. Thus, while the global properties may explain the reactivity, for understanding selectivity the local quantities come into the picture. By direct computation of the local parameters one can probe the sensitivities of different sites in a molecule. The idea is labeled the local hard-soft acid-base (HSAB) principle in analyzing the site selectivity in a molecule.

The local softness describes the response of any particular site of a chemical species (in terms of a change in electron density p(r) to any global change in its chemical potential) and is defined as:

𝑠(𝑟) = 𝑑𝒑 𝑑µ (20)
The local softness condensed to an atom site say k, can be written as:

𝑠(𝑟) = 𝑓(𝑟)𝑆 ( 21
)
Where 𝑓(𝑟) is the Fukui function and 𝑆 is a global softness, which is the integral of 𝑠(𝑟)𝑑𝑟.

Thus, as a result of the relation of 𝑠(𝑟) to the Fukui function 𝑓(𝑟), the local softness is a density-functional concept for characterizing a site and carries the information on site selectivity within a molecule contained in the Fukui function and also the information on relative reactivity from molecule to molecule contained in the global softness. The Fukui function may be thought of as the normalized local softness.

An original direct definition of the local hardness starts from the second functional derivative of the Hohenberg Kohn functional 𝐹[𝑝]. This is the sum of the kinetic energy functional T[p]

and the electron repulsion functional 𝑉 ∞ [𝑝] and is defined for all ground-state densities 𝑝.

The second derivative is the hardness kernel, the two-variable

𝜂(𝑟, 𝑟′) = 𝑑 2 𝐹[𝑝] 𝑑𝑝(𝑟)𝑑𝑝(𝑟′) (22)
The local hardness may then be specified as: 𝜂(𝑟) = 1/𝑁 ∫ 𝑝(𝑟) 𝜂(𝑟, 𝑟 ′ ) 𝑑𝑟′ .

However, several definitions of local hardness have been proposed and compared [START_REF] Chattaraj | Local hardness: a critical account[END_REF][START_REF] Ayers | Local hardness equalization: Exploiting the ambiguity[END_REF] . This local index is then not a local quantity in the sense the local softness is, since it does not integrate to the hardness; consequently, its integral over a given region in a molecule won't necessarily give a regional global hardness.

Local softness and hardness combined with the Fukui functions is thus a basic package to evaluate local reactivity and site selectivity. Consequently, the corresponding QSRP modeling include: regioselectivity of chemical reactions [START_REF] Ayers | Local hardness equalization: Exploiting the ambiguity[END_REF][START_REF] Saha | N-Dependence problem of local hardness parameter[END_REF] , reactivity sequences (intramolecular and intermolecular) of carbonyl compounds toward nucleophilic attack [START_REF] Roy | Local softness and hardness based reactivity descriptors for predicting intra-and intermolecular reactivity sequences: Carbonyl compounds[END_REF] and reactivity of inorganic compounds [START_REF] Chatterjee | Reactivity index scale for interaction of heteroatomic molecules with zeolite framework[END_REF][START_REF] Ghosh | Frontier orbital and density functional study of the variation of the hard-soft behavior of monoborane (BH3) and boron trifluoride (BF3) as a function of angles of reorganization from planar (D-3h) to pyramidal (C-3v) shape[END_REF] .

Electrotopological indices

In contrast to the computationally expensive quantum-chemical descriptors, Kier and Hall 96- 98 provided an easier approach to analyse the molecular structure at the atomic level. The descriptors, called electrotopological state indices and encode the electronic as well as the topological description of the individual constituent atoms of the molecule [START_REF] Hall | Electrotopological state indexes for atom types -a novel combination of electronic, topological, and valence state information[END_REF] , are defined as :

𝑆 = 𝐼 𝑖 + ∆𝐼 𝑖 = 𝐼 𝑖 + ∑ 𝐼 𝑗 -𝐼 𝑖 (𝑑 𝑖𝑗 + 1) 2 𝐴 𝑗=1 ( 23
)
where A is the number of atoms, 𝑑 𝑖𝑗 is the topological distance between the ith and the jth atoms, 𝐼 𝑖 is the intrinsic state of atom, and ∆𝐼 𝑖 is a perturbation factor determined by the influence of the electronic field of a molecule on a particular atom in the molecule. The intrinsic state of the atom is defined as:

𝐼 𝑖 = (2 𝐿𝑖) ⁄ 2 𝛿 𝑖 𝑣 + 1 𝛿 𝑖 ( 24
)
where 𝐿𝑖 is the principal quantum number for atom i, 𝛿 𝑖 𝑣 is the number of valence electrons and 𝛿 𝑖 is the number of sigma electrons. In terms of E-state determination, each atom has its pure intrinsic state perturbed by the electronic environment of every other atom in the molecule. Thus, the intrinsic state encodes the electronic feature of the atom throughout the embodying of the valence electrons which are the most reactive and involved in chemical reactions and bond formations. Further, the presence of the principal quantum number in the expression reflects the differences in the electronegativity of the atoms while the adjacency count of the atom is used to determine its topological features. The ratio of p-and lone-pair electrons over the count of the valence electrons reflects the electronic accessibility and richness of the atom and hence indicated a capability to be involved in intermolecular interactions.

In addition to the individual topological information of an atom, the E-state index can also be utilized to determine the overall contribution of a particular atomic fragment. These indices, called E-state parameters, include the valence state of the atoms of the group along with its hybridization characteristics. Due to the full electronic and topological representation of the group, the E-state parameters are valuable in distinguishing the influence of a certain structural group to the activity profile of the molecule.

Due to its universality and simplicity, the original E-state gave raise to voluminous series of modifications adapted for specific tasks. Thus, the necessity of separated treatment of heavy atoms and their bonded hydrogens for molecules with highly polar functional groups determined the development of hydrogen electrotopological state indices (HE-state indices) complimenting the original E-state indices with electronic and tolopogical information about the chosen hydrogens. Further, this combination has been used for calculation of molecular interaction fields with the assumption that the E-state is defined by superimposing 3D fixed grid over the molecule and hence calculated at each of kth grid point 100 . Generalizing the Estate indices, an additional parameter encoding topological and electronic information related to particular atom types has been proposed as the corresponding sum of E-state values of all atoms of the same atom type in the molecule 101 . Finally, a parameter for bond specification based on the akin bond intrinsic state summed with its perturbation term, has been tentatively proposed for aimed local bond description 102 .

The E-state indices are valuable tools for QSAR studies of biological activities. Thus, it has been successfully applied for the modeling of antithyroid agents with fewer side effects 103 , mutagenicity of aromatic and heteroaromatic amines 104 , anti-inflammatory activity of corticosteroids 105 and receptor binding affinity of progestagens 106 as well as for modeling of fundamental properties such as aqueous solubility 107 and logP coefficient 108 .

ISIDA fragment descriptors

ISIDA descriptors, the development of the Laboratory of Chemoinformatics in Strasbourg, represent specific fragments/substructures present in a molecule [109][110][111][112] and enhanced with the possibility of explicit labeling of an active site. Each substructure is associated to an element i in the descriptor vector, whereas its occurrence in a molecular graph is used as the descriptor value Di. ISIDA fragments could vary in length, topology and inclusion/omission of additional options. The topological variation included sequences of atoms and/or bonds and augmented fragments centered on a certain atom and branching out into concentric circles. The length of the fragments could vary from 1 (the descriptor elements associated to fragments of length 1 are standing for atom counts) to a user-defined number, meaningful for a particular dataset.

As a rule, the size of choice should not be larger than an average molecule's size. The options that could be added to the main description included:

• Formal Charge, permitting to add the information about the formal charge on an atom behind its symbol in the fragment: N +1 -C-N.

• Atom Pairs, types of fragments where two terminal atoms are kept only, with the corresponding topological distance between them: S-6-0, N-3-C.

• All paths exploration, enumerate all the possibilities of the paths between two atoms.

• Dynamic Charge (CGR-specific option), encode local change in charge of the active center atoms while chemical reaction.

The ISIDA descriptors propose two mechanisms that allow the "highlighting" of specific atoms or groups of atoms. The first requires an explicit labeling by the user of the "special" atoms, this is Marked Atom (MA) strategy. The second mechanism exploits the special status of 'dynamic' bonds in Condensed Graph of Reaction (CGR). The CGR-based fragments are generated for a pseudomolecule (CGR), incorporating (condensing) the structures of all the reagents and products. A reaction center is specified by means of special edges, that stand for 'dynamic' (broken/formed) bonds. Figure 3 illustrates the concepts of MA-and CGRapproaches.

The atoms for MA fragments should be labeled with a special flag in a special field of the input file (SDF). Depending on the task, that could be performed by hand, mapping atoms directly in a chemical editor and saving the SDF, or by editing the unlabeled SDF, or with the help of CGR. For the CGR-based fragment, the corresponding input file contains the denotation of the dynamic bonds, constituting the reaction center. The generation of the CGRs could be done by hand or by a special soft, identifying the dynamic bonds by atom mapping.

The preferences of the two types of description should rely on the specificity of the modeling task, i.e. whether the intrinsic nature of the active atoms should be preserved and taken into account, and on the complexity of a chemical transformation. Thus, the task of intermolecular interactions (Figure 3, top) implied an explicit consideration of donor/acceptor (D/A) function of the active atoms, which determines the choice of MA-based descriptors. The descriptors for donors and acceptors in this case could be generated separately and at the end concatenated altogether, always in the same strict order, e.g donor's descriptors-acceptor's descriptors, forming a descriptor vector representing a particular D:A complex. In such manner, the D/A attribution is preserved and will be taken into account by the machine learning algorithm. With regards to chemical reactions, an indisputable advantage of the CGRs is the allowance to encode the structures of all the reactants (reagents and products) altogether with the description of the structural changes. The order of representation of the reactants does not matter in this case. A demonstrative example of the CGRs application is the reactions of cycloaddition (Figure 3, bottom) which involved multiple bond transformations. It should be noted, however, that structural representation of these reactions could be done with the MA-based descriptors as well, by simultaneous labeling of all the atoms of the active center, but the length of the descriptor vector in this case would be too big. In addition, the attribution of the atoms to diene/dienophile (Diels-Alder case) could be lost, if the order of the reactants in the database in not strict. The CGR-based representation thus allows to encode the structures altogether with the active center transformation in a condensed compact form.

For MA-based fragment descriptors an important parameter that could be varied is the degree of 'locality' of the description. An internal mechanism of regularization of the portion of pure local fraction included into the MA-based descriptors is implemented by four marked atom strategies:

▪ -MA0 strategy generates fragments without introduction of the marked atom labels and thus gives a general representation of a structure ▪ -MA1 exclusively generates fragments that starts or ends with the marked atom and hence contains only local fragments ▪ -MA2 exclusively generates fragments that contain the marked atom and is a pure local strategy as well ▪ -MA3 generates all kind of fragments but the marked atom has an explicit label Figure 4 represent the difference between the strategies.

Correspondingly, MA1 is the subset of MA2 descriptor space, which, in turn, the subset of MA3 descriptors altogether with the MA0 nonlocal descriptors.

CGR-based descriptors possess similar, but more restricted option of locality degree regulation: • CGR0 strategy generates all possible fragments

• CGR1 strategy generates only the fragments that contain at least one dynamic bond.

Unlike the CGR-based descriptors, where the fragments are generated for a pseudomolecule, the generation of MA-based descriptors, is performed for each participant of a chemical reaction separately, i.e. for a single molecule. For the case of a molecule possessing more than one marked atom, the two strategies of description are possible:

1. When a local fragment contains more than one labeled atom 2. When local fragments contain one labeled atom and further concatenated with the fragments containing another labeled atom

The preference depends on the specificity of the task. For the case of a bifunctional molecule bearing functional groups (atoms) G1 and G2, the molecule could be represented by local descriptors including both labels (first strategy) or formally represented by two distinct descriptor species -one with focus/label on G1, the other with focus/label on G2. The second approach implicates emphasis on different nature and functions of the groups G1 and G2. For the case of S N 1 dissociation (Figure 5, bottom), the 'active site' is the two atoms connected to the breaking bond, in this case there is no special meaning for these atoms but designating the place of splitting. The fragments will include two labels simultaneously. No difference in labels attribution is encoded. The second instance is hydrogen-bond forming molecule (Figure 5, top) with two binding centers, the first of which is the donor of hydrogen (the corresponding atom is denoted with a green star) and the second is the acceptor of hydrogen atom (denoted with red star). Both atoms are oxygens, thus their nature and their functions during the chemical process should be explicitly designated. These atoms are labeled separately, one by one, and the fragments, generated for each of the atoms then concatenated with each other, so that the atom's nature are encoded by the descriptors number, e.g the descriptors from 1 to 100 reflects the fragments including acceptor atom, the remain descriptors 101 -200 are the fragments with the donor part. In this manner, the atom's functions are preserved and segregated from each other.

ISIDA property-labeled descriptors. Apart from the explicit labeling of active atoms/groups, the ISIDA package supports other, more specific, property labeling [START_REF] Ruggiu | ISIDA Property-Labelled Fragment Descriptors[END_REF]113 . It bears a peculiar information about a particular atom/site that could be maintain alone or coupled with the other fragment generation schemes. These includes:

• partial charge increment

• logP increment

• topological electrostatic potential coloration On the considered example, the histogram of the corresponding property is constructed to estimate the boundaries of the property spectrum to which the atoms are further assigned in accordance to the value of its property, calculated by ChemAxon 114 . Figure 6 (created by MarvinView 115 ), gives an example of logP increment, where the corresponding boundaries and atoms they cover are differentiated by color.

Partial charge labels in ISIDA descriptors calculated according to Gasteiger's method based on the electronegativity of the σand π-bonds. For the logP coloration, the Ghose- Grippen 116 approach is used, according to which the atoms are classified into 120 categories according to their element, oxidation state and the surrounding atoms. The topological electrostatic potential V i on each atom i are calculated from the partial charges according to the formula:

𝑉 𝑖 = 𝑞 𝑖 𝑑 0 ⁄ + ∑ 𝑞 𝑗 𝑑 𝑖𝑗 ⁄ 𝑗≠𝑖
, where 𝑞 𝑖,𝑗 are partial charges on atoms i and j, 𝑑 𝑖𝑗 is the topological distance, and 𝑑 0 is empirically determined virtual distance to take into account the concerned atom charges.

ISIDA fragment descriptors thus propose an efficient way of structural representation. However, it should be noted that relatively long descriptor vector, coming from an exhaustive molecular description, could somehow restricts its applicability.

Descriptors of the reaction conditions

For most of chemical processes, the influence of the experimental conditions is as important as the structural impact of the reagents. Thus, the description of the process should explicitly include the parameters of the solvent medium, temperature, pressure, etc. If for the latter two the description implies just the corresponding magnitude, the solvent description is more complex and needs to take into account the specific and nonspecific solvent effects, able to affect the stability of the formed complex or the transition state, and thus influencing the property.

Accounting for solvent effects is a conventional task in the field of quantum-chemical modeling (QM) and molecular mechanics (MM), where the thermodynamic calculations of chemical, biological or environmental processes referred to the solvated condensed phases.

The solvent models are classified into implicit, explicit and hybrid ones. The earliest attempts of a solvent effect modeling give rise to the implicit models, where the solvent molecules are accounted as isotropic polarizable media. Among widely used, the GB/SA 117 , PCM 118 and COSMO [119][120] approaches should be mentioned. These models are computationally cheap and often provide an acceptable estimation of solvent influence on the process, however, they can not correctly describe specific interactions (e.g. H-bonds), that is of major importance for particular solvents (e.g. water), and thus could fail to describe the influence of media in some systems. More elaborated description is provided by explicit models [121][122][123] , which take the solvent molecules into consideration explicitly, so that, their coordinates and some degrees of freedom are included. These models (AMOEBA, SIBFA, COS) are mostly used in molecular mechanics or molecular dynamics simulation. The parametrization and fitting parameters of these models are therefore derived for a certain solvent group, which, as a consequence, could lead to inability to reproduce some experimental results. The main shortcoming of the explicit methods though is its computational demand. The hybrid QM/MM methods [124][125] , as follows, incorporate the implicit and explicit approaches, so as to provide a reasonable accuracy at fair computational costs. In the frame of the hybrid methods, the energy of the system is composed from the QM-derived energies of the closest to the solute molecular environment, the MM-derived energies of the distant zones and the correction term, refer to the interaction QM/MM energy. The latter is the weakest part of the approach, determining the emergence of various of methods and specific parametrizations for the interaction energy calculation, that complicates the application of the approach. The hybrid methods thus could be tested and used thoroughly, so as to give a reliable result.

In QSPR, it is important to provide the machine learner with enough information about the solvent properties. Machine learning will then figure out which of the provided solvent properties appear to correlate with the modeled property of the studied solute(s), and propose the optimal (non)linear functional form to express the dependence of the latter with respect to the former. So far, however, there were too few works related to the modeling of chemical processes with the account for the solvent environment. The one that explicitly include the solvent parameters is related to the modeling of S N 2 reaction 126 . The solvent media has been rendered there with the six parameters, standing for polarity and polarizability, derived on the basis of the works of Born 127 and Kirkwood 128 . These works emphasize the importance of nonspecific solvation, which is determined by polarity and polarizability. The former term could be expressed by the three functions of dielectric constant ɛ 129 , whereas the polarizability could be represented by the three functions of refractive index n D 20 129 .

The consideration of the solvent effects, however, should be completed with the inclusion of the specific solvation term. These parameters should meet the following requirement: derived from the experiment, be measured for a large set of different solvents and not to be biased by the probe, used during the experiment, but referred to a 'general' solute behavior toward a particular solvent. These preconditions are met for the solvents effect scales, among which Kamlet-Taft solvent effect scale is one of the most widely known. The first one, α scale 130 is referred to hydrogen bonding ability of solvents. This scale is based on solvatochromic parameters, averaged for several probes, so that it has a built-in 'fuzziness' and measure the ability to donate hydrogen bonds of the solvent molecules to a 'general solute', rather than specifically for the probe employed in the experiment. The second one,  scale 131 , is based on the ultraviolet-visible spectral band of suitable probes. This is again an averaged quantity, for which the wavenumber shifts of several protic indicators, relative to structurally similar, but aprotic probes, are used. The third value is the π* value 132 , based on the average of values of the π →π* transition energies for several nitro-substituted aromatic indicators. The quantity is normalized to give pi* = 0 to cyclohexane and pi* = 1 to dimethylsulfoxide, and, as for the previous two scales, multiple probes are used to eliminate specific interactions and spectral anomalies. This value measures a certain 'blend' of polarity and polarizability. Another set of specific interaction-referred scale is the Catalan parameters. Similar to Kamlet-Taft scales, the Catalan parameters included the solvent polarity/ polarizability (SPP) 133 , solvent basicity (SB) 133 and solvent acidity (SA) 134 measures.

In the frame of this thesis, the solvent effects have been described with the following set of the experimental parameters:

• four functions of dielectric constant ε, standing for nonspecific interactions : Born

  1 - = B f , Kirkwood 1 2 1 + - =   K f , 1 1 1 + - =   f and 2 1 2 + - =   f .
• three functions of the refractive index nD [START_REF] Roberts | σ-constants of the carbethoxyl and hydroxyl groups[END_REF] (n), as well reflecting the nonspecific effects:

2 1 2 2 1 + - = n n g , 1 2
1 2 2 2 + - = n n g , ) 1 2 )( 1 2 ( ) 1 )( 1 ( 2 2 + + - - =   n n h .
• Kamlet-Taft's α, β and π* parameters • Catalan's SPP, SA and SB constants

An additional reaction condition term of temperature has been included as the reversed value:

1/T in Kelvin degrees, as it corresponds to the Arrhenius equation (ln 𝑘 = ln 𝐴 -𝐸𝑎 𝑅𝑇 ⁄ ). In case of water-organic mixtures, a molar fraction of the organic solvent has been added as a descriptor. The hierarchical clustering dendrogram of the solvents used in different projects of this study and based on the 13 chosen solvent parameters is given on Figure 7. QSPR modeling is the determination of a mathematical relationship between the chemical structure (or more complex information included in the description) and a modeled property/activity:

𝑷𝒓𝒐𝒑𝒆𝒓𝒕𝒚 (𝑎𝑐𝑡𝑖𝑣𝑖𝑡𝑦) = 𝒇(𝑐ℎ𝑒𝑚𝑖𝑐𝑎𝑙 𝑠𝑡𝑟𝑢𝑐𝑡𝑢𝑟𝑒, 𝑠𝑜𝑙𝑣𝑒𝑛𝑡, 𝑇, 𝑝𝑟𝑒𝑠𝑠𝑢𝑟𝑒, 𝑒𝑡𝑠)
The corresponding structure-related information in the function argument should be encoded numerically, composing the descriptor vector, defining the position of an object in chemical space. The object's numerical description then used to build a QSPR model, embodied the mentioned mathematical relationship, the goal of which is to be able to predict a certain property/activity over a wide range of new (in a sense they were not used for a model's building) chemical objects.

A general used for model building procedure of QSPR model building included the following steps:

i. Data curation and standardizations. This includes rejection of entries associated to missing or chemically invalid structures, missing or unreliable experimental endpoints, removal of counterions not needed for modeling, conversion of structures into a common, standardized representation style (aromatic bonds, split-charge nitro groups, etc) and removal of duplicates (complete removal if a same structure is associated to conflictingly different experimental values). ii. Descriptors calculation. In this work, we used ISIDA fragment (Marked Atom or CGR-based) descriptors for structural representation of a chemical object (described in section 2. 1.4). Accordingly, the descriptor vector is constructed from structural fragments of different length with the corresponding occurrences being the descriptor's value. The descriptors' values are then normalized form 0 to 1. iii. Model building. The algorithms used for building of the predictive models called the machine learning methods. Whereas there is a diversity of different algorithms, the ones chosen for modeling in this work are the state-of-the-art techniques proved its efficiency and usability: the Support Vector Machines and Multiple Linear Regression. The tool combining both possibilities, regression and visualization, is the Generative Topographic Mapping. The description of these methods is given in section 3.2. A data set used for building of the model is called training set iv. Model validation. This includes cross-validation, performed at the stage of model building (cross-validation is described in section 3.3.1) and evaluating the model performance on the training set and external validation. For the latter, some part of data could be excluded from the initial data set at the beginning or new data could be used. These data were not involved in model's building at all. This test set is required to assess predictive power of the model and its utility.

The constructed model is used for the prediction of the corresponding property/activity of unknown or untested structures thus providing, of course depending on the overall accuracy of the model, a numeric value or a classification belonging of the structures given.

Machine Learning algorithms

Support Vector Machine (SVM)

Firstly introduced by V.Vapnik [135][136] , Support Vector Machine at its origin is a binary classifier that finds a separating hyperplane so as to best segregate the two classes (Figure 8). The hyperplane is constructed in such a way to have the biggest gap between the hyperplane and the data instances of either side, so that to minimize the chance of misclassification.

Accordingly, the SVM model then categorizes new examples according to which side of the hyperplane they fall.

The method however is suitable not only for linearly separated objects, but for nonlinear task as well. To perform the separation in this case, a new coordinate is introduced in such a way that in the resulting higher-dimensional space the classes are easily separated (Figure 9). The possibility to perform nonlinear separation in case of usage of linear classifier utilizes so called kernel trick and is widely used in kernel-based methods. The Support Vector Machines for the regression task implies optimization of the regression function f(x), which is searched under the following constraint:

|𝑦 -𝑓(𝑥)| < 𝜀 ( 25
)
here 𝜀 is the error threshold. The function is optimizing till no errors largen than 𝜀 are produced. However, to be able to model data with persisting errors, a proportional cost is introduced:

𝜉 = |𝑦 -𝑓(𝑥) -𝜀| ( 26
)
The model is fitted so that to minimize its complexity and the proportional cost.

In this work the Support Vector Regression (SVR) is the main or benchmarking method of modeling, the developed models of which are the resulting web-implemented output of the projects.

Multiple Linear Regression (MLR)

Multiple Linear Regression aimed at finding the equation between the property/activity and the descriptors, encoded the chemical object, with a crucial assumption that the relationship is linear. That could be represented as follows:

𝑷𝒓𝒐𝒑𝒆𝒓𝒕𝒚 = 𝛼 + 𝛽 1 * 𝑑𝑒𝑠𝑐𝑟𝑖𝑝𝑡𝑜𝑟1 + 𝛽 2 * 𝑑𝑒𝑠𝑐𝑟𝑖𝑝𝑡𝑜𝑟2 … + 𝛽 𝑁 * 𝑑𝑒𝑠𝑐𝑟𝑖𝑝𝑡𝑜𝑟𝑁
The goal thus is to find the intercept 𝛼 and the corresponding regression coefficient 𝛽 𝑥 , which could be considered as a contribution of a certain descriptor into the property i.e. measures the unit change in the dependent variable with the change of the descriptor, so as to fit the property variation.

The software for MLR in this study is ISIDA QSPR 1 , which combines backward and forward stepwise variable selection (prior selection of those variables influencing most on the model's predictive ability, needed to provide more robust and cost-effective prediction) generating a large number of linear models forward by the selection of the most robust ones for the consensus prediction, that is an average of the estimated property values obtained with the selected individual models.

Generative Topographic Mapping (GTM)

Generative Topographic Mapping is a method combining the modeling capability along with data visualization and data analysis tools. Firstly introduced by Bishop 3 in 1990th, the method performs a non-linear dimensionality reduction of the D-dimensional data space (where D is the number of descriptors) onto a 2-dimesional latent space (GTM map) by embedding a flexible 2D manifold into the D-dimensional data (Figure 10).

The 'unfolded' manifold is a square grid of K nodes. The assignment (mapping) of the nodes to the manifold points is defined by a mapping function 𝑦 𝑑 (𝑥; 𝑊) set up with the help of M radial basis functions (RBFs).

𝑦 𝑑 (𝑥; 𝑊) = ∑ 𝑊 𝑚𝑑 exp ( ‖𝑥 -𝑥 𝑚 ‖ 2 2𝜎 ) 𝑀 𝑚=1 (27) where d goes from 1 to D, W is the 𝑀 × 𝐷 weight matrix connecting RBF and data space points, 𝑥 𝑚 is the center of the m-th RBF. The overall number of the RBFs and the width σ are the optimizable parameters of the method.

As a probabilistic extension of the Kohonen Self-Organizing Maps, GTM operates with probabilities of a data object to be mapped into a certain node. Moreover, the object has nonzero probability over all the nodes. Consequently, it could be characterized by its probabilities, which are called responsibilities. Responsibilities constitute a responsibility vector, the main descriptive characteristic of an object, used for class belonging assignation in case of classification, or property value calculation for the case of regression modeling.

The responsibility of the k-th node for n-th data point 𝑡 𝑛 is calculated using Bayes' theorem: 𝐾 𝑘=1 (28) Responsibilities are normalized over the grid of nodes and their sum for a given item is equal to 1.

GTM-based regression and classification models

The general procedure of the model building is similar for regression and classification tasks, as in both cases an object (for simple case -a single molecule) is characterized by the assigned responsibility vector. The following steps are included:

i. Obtaining a GTM grid, each node of which will be attributed with the corresponding responsibilities of every molecule of the training set. ii. Defining the property/class value of each node based on the contribution of each molecule to a certain node (which is molecule's responsibility) and its corresponding property/class value. This procedure is called coloration, since the property distribution will be expressed as the color profile of the map. The same GTM map hence could be colored differently depending on the training set property/class values. iii. Projecting the test set compounds and calculating its responsibility vectors further used for the property/class prediction.

The node property value 𝐴 ̂𝑘 of GTM regression model (step ii) is calculated as follow:

𝐴 ̂𝑘 = ∑ 𝐴 𝑛 𝑅 𝑘𝑛 𝑁 𝑛=1 ∑ 𝑅 𝑘𝑛 𝑁 𝑛=1 ( 29
)
where N is the number of molecules, 𝐴 𝑛 is the experimental property of the n-th molecule, 𝑅 𝑘𝑛 is its responsibility in the k-th node (see eq. 28). The calculated node property values are used for the GTM activity landscape representation (eq. 29.

2) -the final result of the map coloration (detailed below).

𝐴 𝑞 = ∑ 𝐴 ̂𝑘𝑅 𝑘𝑞 (29.2)
Similarly, GTM classification model attributes a class accessory to each node by averaging the responsibilities over all training set compounds, then, for any q-th test set compound, the probability to belong to the i-th class 𝑃(𝐶 𝑖 |𝑞) is calculated according to the formulae:

𝑃(𝐶 𝑖 |𝑞) = ∑ 𝑃(𝐶 𝑖 |𝑘) × 𝑅 𝑞𝑘 𝑘 ( 30
)
where P(C i |k) is the conditional probability of the class C i for the given node k, calculated according to the Bayes' theorem:

𝑃(𝐶 𝑖 |𝑘) = 𝑃(𝑘|𝐶 𝑖 ) × 𝑃(𝐶 𝑖 ) ∑ 𝑃(𝑘|𝐶 𝑗 ) × 𝑃(𝐶 𝑗 )
𝐶 𝑗 (31) where 𝑃(𝐶 𝑖 ) and 𝑃(𝑘|𝐶 𝑖 ) are, respectively, a fraction of compounds of the class C i and a normalized cumulated responsibility of the class C i in the training set.

GTM visualization

The ability to visualize the data distribution is the main advantage of GTM over the classical machine learning methods. The method proposes different schemes of representation of data distribution thus allowing to analyze various aspects of data. The main visualization techniques used herein are GTM property landscape and GTM class landscape, correspondingly, representing data form regression or classification side. The landscapes are created according to the mentioned equations (eq. 29-31) and reflect the node's property-or class attribution. In addition, the landscapes are weighted by data density: the more molecules (or more complex chemical objects) are located near a certain node, the more opaque the color of the node, correspondingly, if no molecules are projected into a node with any reasonable responsibility, the node remains to be transparent (blank). Figure 11 shows examples of property-and class landscapes providing the equilibrium constant distribution and the tautomeric type separation of 695 tautomeric equilibria (the project is described in section chapter 6).

Model quality estimation

Cross-validation and external validation

A QSPR model needs to be validated in order to estimate whether it has a key competencethe ability to predict the property of new objects. The model performance should not be evaluated on the data that was used to build the model: indeed, the model would just repeat the property value/class of the samples that it has just seen and shows a perfect score (if the samples were not 'outliers' constantly mispredicted by the model), but would fail to predict yet-unseen data. This situation is called overfitting. To avoid overfitting, the model's performance is estimated during cross-validation, the procedure that envisages retention of a part of data and its further usage for model's evaluation. The initial data is thus divided into two parts, training set and test set. However, to get unbiased independent predictions for each object of the data, all of them should be estimated during test set prediction. To do so, the portioning of data is performed several times, usually 5, times, correspondingly called 5fold cross-validation. Each time the different 5 th part of data is retained as a test set, and the other 4/5 are used as the training set (Figure 12). That insures the unbiased prediction will be obtained for each object of the data.

As a rule, the built model should be also estimated on a data set not at all related to the initial data (used for building) and comes from different source, or, if not possible, randomly chosen from the initial data before any modeling, and retained. This set, called external test set, is a key tool for model's performance analysis and shortcomings revealing.

Regression-and classification model's performance criteria

The predictive performance of a regression model (estimating continuous property) is obtained with the following parameters:

▪ Determination coefficient:

𝑄 2 (𝑜𝑟 𝑅 2 ) = 1 -∑(𝑌 𝑒𝑥𝑝,𝑖 -𝑌 𝑝𝑟𝑒𝑑,𝑖 ) 2 / ∑(𝑌 𝑒𝑥𝑝,𝑖 -< 𝑌 > 𝑒𝑥𝑝 ) 2 𝑛 𝑖=1 𝑛 𝑖=1 (32)
▪ Root Mean Squared Error: where Yexp and Ypred are, respectively experimental and predicted values of property and <Y>exp is the mean of experimental values.

𝑅𝑀𝑆𝐸 = [1/𝑛 ∑(𝑌 𝑒𝑥𝑝,𝑖 -𝑌 𝑝𝑟𝑒𝑑,𝑖 ) 2 𝑛 𝑖=1 ] 1/2 (33)
Classification models (predicting the label of an object, i.e. active/inactive) estimated here by the following:

▪ True Positive Rate:

𝑇𝑃𝑅 = 𝑇𝑃/𝑃 ( 34
)
where TP is the number of True Positive (being positive and predicted as positive) species while P is the overall number of experimentally positive class species in the data set. ▪ True Negative Rate:

𝑇𝑁𝑅 = 𝑇𝑁/𝑁 ( 35
)
where, similarly, TN is the number of True Negative (being negative and predicted as negative as well) species and N is the number of experimentally negative class species. ▪ Balanced Accuracy:

𝐵𝐴 = 𝑇𝑃𝑅 + 𝑇𝑁𝑅 2 (36)

Applicability Domain

Applicability Domain (AD) defines the area of chemical space where the model is presumably accurate. The concept of Applicability Domain assumes that the objects similar to those used for model building, will be predicted accurately rather than very different, in terms of descriptor vector similarity, targets. Figure 13 gives an illustration on the example of 2D chemical space.

PART II. RESULTS AND DISCUSSIONS

Chapter 4

QSPR modeling of halogen bond basicity of binding sites of polyfunctional molecules

Halogen bonding started to attract specific attention across the chemical, biochemical, and material sciences very recently -this peculiar interaction received its official IUPAC definition only in 2013 137 . Indeed, a predisposition of halogen atoms to behave as nucleophiles due to their high electronegativity is a well-established understanding, thereupon halogen atoms are mostly considered as the regions of high electron density. However, their ability to behave as Lewis acids notwithstanding their intrinsic nature was revealed in the beginning of the 1900 th by the formation of complexes such as Hal 2 …NH 3 and Hal 2 …OH 2 138-139 . Nonetheless, the revelation of the electron density in halogen atoms being anisotropically distributed whenever the atom is covalently bound to one or more atoms, has emerged only recently [140][141] .

A covalently-bonded halogen atom surrounded by the area of rich electron density forming a belt, orthogonal to the covalent bond, where the electrostatic potential is negative, but at the same time the electronic distribution anisotropy shapes a region of lower electron density (the so-called σ-hole) where the potential is frequently positive. This region can form attractive interactions with electron-rich sites, determining the ability of halogen atoms to interact with nucleophiles. Figure 14 gives an illustration, where the color code corresponds to the value of the electrostatic potential at the surface. Accordingly, the scale of halogen bond atoms' strength to act as Lewis acids is referred as follows: 𝐹 < 𝐶𝑙 < 𝐵𝑟 < 𝐼. Fluorine, as less polarizable one, is less prone to participate in halogen bonding and being capable of one only when attached to particularly strong electronwithdrawing groups. Iodine therefore is the most active and convenient for experimental studies.

The role of halogen bonds (XB) is particularly prominent in the areas of crystal engineering 142 , but also important for the elaboration of three-dimensional networks and the formation of liquid crystal phases, in different areas such as biological molecules design [143][144] and nanotechnologies 140 . The comprehensive outline of recent advances and historical perspective of the field are reviewed in works of Cavallo 145 , works of Legon 146 and Priimagi 147 .

Our aim in this project was the development of the universal scale of halogen bond acceptor strength, i.e. halogen bond basicity, that could be considered as a scale of nucleophilicity as well. The basicity scale based on the strength of the complexation with diiodine, was the object of our publication in Molecular Informatics 148 . The asset of the paper is the efficiency of the developed scale for the prediction of halogen bond strength of not only monofunctional, but polyfunctional species as well, that expands its applicability toward complex biological molecules and supramolecular building blocks. Due to their low computational costs, the developed models are of practical relevance for an efficient screening of large sets of compounds. In the paper we also discuss the borderline of the applicability of the constructed scale, providing the examples of molecular species (article's section 3.2) possessing structural features responsible for certain steric effects, affecting the complexation constant, that should be thus treated cautiously. Apart from the contribution toward halogen-bond driven processes, one could find the section of comparison of the strength of halogen and hydrogen bonding to be of particular interest (article's section 3.4).

The project was done in collaboration with Jerome Graton and Jean-Yves Le Questel (University of Nantes) providing thorough experimental data, and Vitaly Solov'ev (Institute of Physical Chemistry, Moscow) conducting the MLR-related part of calculation as well as the effective constant evaluation by means of ChemEqui 149 software.

The corresponding article is given in the end of the section, with the authorization of all authors.

Modeled object and property

Herein, the halogen bond (XB) donor molecule is the same for all the complexes, hence, its structure could be excluded from consideration. The modeled object is thus the structure of a designated XB acceptor, the active center of which, binding the halogen atom, is attributed.

The modeled property is the complexation constant (log KB I2 ) of an organic molecule, considered as a Lewis base, with diiodine (I 2 ). The experimental values are referred to 1:1 complexation in hexane at 298K. The structure of a XB acceptor is represented by the Marked Atom (MA)-based descriptors, where the corresponding marked atom is attributed to the active site of a molecule, binding with I 2. All four marked atom strategies of MA-based descriptors have been tried and compared. In case if molecules possessing two putative binding centers, the main one has been indicated in the initial source.

Modeling workflow 4.3 Data preparation

Initial data has been reported in the work of Laurence et al. 150 . The log KB I2 values represented primary the experimentally measured values in heptane, hexane, cyclohexane and methylcyclohexane at 298K. Their differences as a result of solvent effects are generally within the experimental uncertainties 151 . Some of the log KB I2 have been experimentally measured in CCl4, CH2Cl2 or CHCl3. In this case, the corresponding values have been recalculated and referred to hexane by known linear Gibbs energy relationship 151 . If a compound had several reported equilibrium constants values in different solvents, only the primary value in alkane was selected. Cis/trans-isomers with diverging XB acceptor propensities have been removed. The structures have been standardized according to the procedure used on our virtual screening web server (http://infochim.u-strasbg.fr/webserv) and based on ChemAxon's Standardizer 152 (neutralization, isotopes removal, conversion to 'basic' aromatic form etc.) The labeling of the concerned XB active sites, explicitly assigned in the initial database, has been performed manually. Thus, the training set consisted of 598 organic molecules of 14 different types of XB acceptor atoms, the weakest of which is the π-electronic carbon and the strongest is the sulfur of thiophosphoryl group (App., part I, Table I.1).

Figure 15. Workflow of the modeling of the strength of halogen bonding between organic accceptors and diiodine (logKBI2).

The external set consisted of 11 polyfunctional species, collected from the same source, for which the effective constants (log K B eff), attributed to the binding that involved all of the active centers, have been experimentally measured. The measurements performed in solvents, different from hexane, have been recalculated to the appropriate solvent by known 151 linear free energy relationships.

Computational details

The fragmentation schemes included various atom coloration by elements symbol, by CVFF force field label or by pharmacophoric types (described in section 2. 1.4). The considered fragment topologies were sequences and atom-centered fragments of the minimal length from 2 to 4 and the maximal length from 3 to 8. Overall, 480 descriptor sets (120 for each marked atom strategy) have been tested. The modeling has been performed by SVR and MLR methods. The performance of the models has been estimated by the R 2 and RMSE values in 5-fold cross-validation. The applicability domain control method was Bounding Box. The most robust SVR and MLR models constituted the consensus SVR (MLR) models (CM), rendering the property value as the corresponding average of the values, predicted by the individual models. The prepared consensus SVR and MLR models have been further used for the prediction of the external test set.

The assessment of log K B eff values for polyfunctional molecules of the external set was derived from the predicted log K BI2 of each individual active center of a molecule. The corresponding estimation was done with the help of ChemEqui 153 program simulating a network of chemical equilibria in solution and designed to handle the cases of simultaneous coexistence of several mono-and polybinded species.

The details of the computational procedure including the specification of the scanned descriptor spaces for SVR and MLR, list of the models included into the final web-deployed consensus SVR model, as well as the detailed workflow of data curation and treatment are described in the article.

Results and discussions 4.4.1 Cross-validation.

The average statistical values in 5-fold cross validation returned by the consensus SVR/MLR prediction for each marked atom strategy are summarized in Table 2. The best result is achieved by the MA3 strategy, which explicitly distinguish fragments belonging to the reaction center and its environment. The MA3 strategy could be seen as the sum of MA2 and MA0 descriptors, where the MA2 describes the immediate surrounding of the active center and MA0 does not pinpoint the active center but generally characterizes the molecule. The MA3 strategy thus has been chosen for the construction of the consensus SVR and MLR models (the constituting individual models are listed in Appendix, part I, Table I.2) which were used for the prediction of the external set of 11 polyfunctional molecules.

External validation.

The predicted values obtained by the model for individual binding centers were combined into the effective constants (log K B eff), the comparison of which with the experimental values is given on the Figure 16. The predicted effective constant reproduces the experimental logK B eff with the RMSE values close to the ones of the cross-validation stage.

Comparison of the strength of halogen and hydrogen bonding

A set of 166 molecules with available H-bonding acceptor strength data against 4-F-phenol in CCl 4 154 has been used for the evaluation of the relationship between halogen and hydrogen bonding strength. As expected, the general overall tendency, concerning the predominance of one or another binding type, could not be observed, however, strong correlations have been found within specific chemical families (Table 3; App., part I, Table I This can be explained by the fact that the physico-chemical nature of the acceptor sitechemical element, hybridization, etc -defines the generic order of the magnitude of the interactions. Within a given family, the generic acceptor propensity of the center is modulated by its chemical context -and it turns out that this modulating effect is comparable for both H-bonding and XB strength: a same substituent will impact both properties similarly. As it was determined, for oxygen bases (the C=O, -O-, P=O, S=O sites), the log K BI2 is lower than the log K BHX of the H-bond complexes. An opposite trend is observed for sulfur bases (the -S-, C=S, P=S sites), for which halogen-bonding is considerably stronger than Hbonding. Similar regularity is observed for the stability of diiodine complexes with primary, secondary and tertiary amines, as well as in the case of aromatic nitrogen bases, for which the stability of diiodine complexation is compatible or higher than the stability of H-bond complexation.

Conclusion

This project is a starting point of this thesis, representing the simplest case of the modeling for the system, where one of the participants and the experimental conditions stay constant. The modeled property is the binding strength of complexes of organic molecules with diiodine (I 2 ). The quantitative value of the binding strength serves as the halogen bond basicity scale, or, more general, as a scale of nucleophilicity of organic molecules. Here we report a successful QSPR modeling of the halogen bond basicity of 598 organic molecules for which the binding constants have been measured at unified conditions (hexane, 298K). The structure of an organic molecule has been characterized by Marked Atom-based descriptors, the different labeling strategies of which, representing particular levels of generalization/specification of structural description have been applied and compared. The MA3 strategy turned out to be the best performer as it combines an explicit characterization of the active sites with the description of the overall structural arrangement of a molecule. The cross-validation results of the SVR and MLR individual models built on the MA3-based descriptors spaces are close to the experimental errors: RMSE=0.39-0.43 (R 2 =0.92-0.93).

That should be noted, that regardless of the use of the best fragment descriptors, these can not cover the entire range of different structural and electronic effects playing a role in the complexation strength, notably the bidentate halogen/hydrogen bond interaction scenarios occurring in certain conformations. However, during model building, the fitting errors caused by any types of similar effects are minimized so that their average affect over the property is below the intrinsic imprecision of the model. Extensively cross-validated consensus SVR and MLR models have been challenged to predict the effective complexation constants of polyfunctional molecules of the external test set. The models showed robust cross-validation statistics (R2 =0.70-0.81, RMSE=0.44-0.56) and were able to successfully extrapolate the interaction of polyfunctional compounds with I 2 , for which the experimental effective binding constant could be inferred from the individual propensities of all the groups putatively participating in XB. The comparison of the H-bond and XB complexation constants does not show any global relationship between these related, but mechanistically quite different chemical interactions. However, strong piecewise correlations within chemical families based on the same type of H-bond/XB donor were found, which means that while the intrinsic HB or XB strength of these centers are uncorrelated, the modulating impact of the substituents on both HB and XB are comparable.

A predictor of the halogen-bond basicity of acceptor sites of organic molecules was created on the entire training set and comprises the best performing SVR models (App., part I, Table I.2). The consensus model is publicly available on the web server: http://infochim.ustrasbg.fr/webserv/VSEngine.html, altogether with the automatic binding centers labeling and molecule's applicability domain estimation.

Chapter 5

QSPR modeling of the Free Energy of hydrogenbonded complexes with single and cooperative hydrogen bonds.

Discovered first around 100 years ago 155 , hydrogen bond is still object of numerous research and debates. The reason of this long-lasting interest is determined by the importance of hydrogen-bond based interactions to a broad spectrum of fields varying from biology to material science. The topic of hydrogen bond interactions drew particular attention in 1990th with the boom in developing of supramolecular and crystal engineering researches. Since that time, the depth and the complexity of the phenomena have expanded drastically. A new concept of hydrogen bonds has been emerged 156 and new aspects of weak hydrogen bonding occurring in biological systems have been discovered 157 . Hydrogen bond in the present time is interpreted as a region alternating from covalent bonds to van der Waals interactions, ionic interaction and even π-cation interchange.

The following definition is proposed for the complexity of hydrogen-based interactions: An X-H … A interaction is called a "hydrogen bond" (HB), if (1) it constitutes a local bond, and (

X-H acts as proton donor to A. The energy range of hydrogen bond dissociation varies from 4 to 160 kJ/mol, and the distance up to 3.2Å is considered potentially capable of bonding. Within this range, the nature of the interaction is not constant, but includes electrostatic, covalent, and dispersion contributions in varying weights. More about the nature and the variety of HB interactions could be found in numerous related books 156,[158][159] , works of Desiraju 160-161 , Leiserowitz 162-163 and Steiner 164 , while the biological aspects could be found in a book of Jeffrey and Saenger 165 .

Following the line of mono-and polyfunctional intermolecular interactions, this project deals with the modeling of the strength of hydrogen-bonded complexes. Our paper on this topic has been published in Molecular Informatics 166 . The task is aimed at the modeling of the strength of both, mono-and polyfunctional hydrogen bonds, which, this time is formed by different acceptors and different donors. The data set used for modeling is so-far the biggest used for QSPR study of hydrogen bond strength. Initial data included the measurements of the same complexes in different solvents, that allows to build linear correlations, so as to go from the reference solvent to the required one, allowing the comparison of HB strength relative to different media. The obtained linear correlations have been estimated during the external validation. The performance of the developed models has been evaluated on two external sets. The first one was formed by the complexes with single HB, among which there where donors/acceptors encountered in the training set as well as structurally unknown molecules. This fact suggests a different from the traditional, ubiquitously applicable for single molecules, manner of the model's predictive performance estimation. Here, the external set complexes have been attributed to four different classes that correspond to a certain degree of 'novelty' of the complex with respect to the training set. Consequently, the model's performance has been evaluated for each of the classes.

The work was done in collaboration with Vitaly Solov'ev (Institute of Physical Chemistry, Moscow) carried out the MLR calculation.

The article is given in the end of this section, with the authorization of all authors.

Modeled object and property

The data set, referred to 'standard' conditions, was composed of varying donors and acceptors coupled by single hydrogen bonds. The active site of each participating molecule was attributed and explicitly marked: correspondingly, this is the donor atom, providing hydrogen, and the acceptor atom with free electron pair. The Marked Atom (MA)-based descriptors have been used. The structures of the donor and the acceptor molecules have been treated separately: for each of the participants the fragment descriptors, including the local ones, i.e. with the corresponding donor/acceptor (D/A) labels on the active sites, have been prepared and in the end concatenated altogether, forming an integrated descriptor vector, representing a particular donor-acceptor complex. The MA3 strategy was the best one in the previous project with similar task of modeling of the strength of intermolecular complexation, which could be explained by the fact that the strategy encompasses and describes the whole molecule but at the same time explicitly distinguish fragments belonging to the reaction center. Thus, the MA3 strategy is involved herein.

Modeling workflow

Data preparation

An initial data set consisted of 4002 HB complexes measured in different solvents at several temperatures has been compiled from the literature 151,154,[167][168][169][170] . The complexes have been

Computational details

Out of four available marked atom strategies, the MA3 strategy was chosen for the modeling, based on a preliminary MLR study (App., part II, Table II Eventually, sets featuring new donors in interactions with new acceptors are still a bigger challenge. Thus, the monofunctional external set could be considered as containing four classes of complexes attributed to four distinct degrees of 'novelty':

• PAIROUT-both donor and acceptor were featured in some of the training set complexes, but never together;

• ACCOUT-the acceptor of this pair was not included the training set, but the donor was present in some HB complexes;

• DONOUT-the donor of this pair was not in the training set, but the acceptor was seen in some HB complexes;

• BOTHOUT-neither donor nor acceptor were in the training set.

The model's predictive performance thus is referred to these four specific classes represent an increasing degree of difficulty of extrapolation. The ΔG assessment for the second external test set with cooperative HBs was calculated with the assumption that the observed experimental affinity linearly correlates with the sum of individually assessed ΔG values.

The performance of the models has been estimated by R 2 and RMSE values. The applicability domain was defined by Bounding Box. Winning SVR and MLR models constituted the consensus SVR/MLR models, that predict the property as an average of the values, predicted by individual SVR (MLR) models.

The details of the computational procedure concerning the specification of the winning descriptor spaces, the list of the descriptor spaces and the corresponding SVR parameters of the individual models included into the final consensus SVR model (CM), as well as the detailed protocol of data curation and treatment are described in the corresponding sections of the article.

Results and discussions

Cross-validation

The performance of the individual models has been estimated in three-fold cross-validation (3CV). Since various descriptors represent various, complementary points of view of the molecular structure, the individual models using them capture the chemical information of different nature. Therefore, the combination of the 40 individual SVR models into the SVR CM model leads to a significant synergetic effect, boosting RMSE to 1.50 and R 2 to 0.94. The performance achieved with consensus MLR calculations (RMSE = 2.22 and R 2 = 0.87), see Figure 19, is less impressive. There may be several reasons for this:

a. The ISIDA MLR tool automatically scans through possible fragmentation schemes, but has no access to the "colored" ISIDA descriptors (described in section 2. 1.4) that were manually added to the pool of candidate SVR descriptor spaces, and were found to win the competition. b. Non-linearity seems to play an important role in HB affinity modeling: albeit the linear kernel was available amongst SVR options, only two models out of the 40 constituting the CM incorporated this option, and both of their ranks are at the bottom of fitness-ranked list. Accounting for model applicability domain slightly improves predictive performance because of discarding some 10% compounds: RMSE MLR (within AD) = 2.11 kJ/mol.

External validation

External test set 1 (complexes with single hydrogen bond)

Results given in Table 4 show that the predictive performance of both SVR and MLR consensus models is not as good as the one observed in cross-validation (Figure 19). This can be explained by the noise caused by the inclusion of solvent corrections as well as the fact that one third of the compounds are outside of the models' applicability domain. Discarding these species resulted in significant decrease of RMSE till 2.5-3.01 kJ/mol. The comparison of four validation scenarios corresponding to different degrees of novelty revealed that the accuracy of the prediction decreases in order: DONOUT > ACCOUT > BOTHOUT > PAIROUT.

External test set 2 (complexes with cooperative hydrogen bonds)

The predicted ΔG values for 12 complexes with cooperative HBs have been assessed from the sum of the Free energies for individual centers (App., part II, Table II.3) according to the formulae:

𝛥𝐺𝑝𝑟𝑒𝑑 = 𝛼 ∑ ∆𝐺 𝑝𝑟𝑒𝑑, 𝑖 𝑛 𝑖=1 ( 37
)
where 𝑛 is the number of HBs (in our case 𝑛 = 2). The parameter 𝛼 has been fitted by the least squares method, and equal to 0.60±0.02 for SVR and 0.66±0.04 for MLR predictions.

The corresponding graphic of the predicted vs experimental values is given in Figure 20. Thus, a reasonable correlation observed for the values with the performance similar to that on the cross-validation stage: RMSE = 1.63 and 2.68 kJ/mol, R 2 = 0.87 and 0.65 for SVR and MLR, respectively.

Conclusion

By contrast to the previous project where one of the interacting entities remained constant, this project is devoted to the modeling of the strength of intermolecular complexes formed by varying donors and acceptors. The structures of the complexes have been represented by the Marked Atom (MA)-based descriptors, where the corresponding atoms that have been labeled are the donor and the acceptor of hydrogen. The strength of hydrogen bond complexation was characterized by the experimental Free Energies (ΔG, kJ/mol) measured at 'standard' conditions: CCl 4 , 298K. To our knowledge, the data set utilized (3373 complexes) is so far the largest used for hydrogen-bond complexation propensity prediction. The cross-validation performances of the models are similar (R 2 =0.87-0. accumulation of errors on the individually predicted pKa terms. Furthermore, the approach does not consider varying reaction conditions. Here we propose to treat the problem of tautomeric equilibria evaluation directly, based on the experimental K T values and Marked Atom-based descriptors characterizing the structure of the molecules and the character of the transformation. The data set included 10 different types of tautomeric transformations measured in different solvents and temperatures. Presence of reaction conditions allows a fullblown modeling of the combined structural and reaction condition impact on the equilibrium constant as well as to seek for the specific patterns characterizing various types of tautomerism.

This is the first time the Generative Topographic Mapping (described in section 3.2. 3) is challenged in modeling and visualization of combined structure/conditions chemical data. An impact of the reaction conditions on the property could be examined by means of GTM maps exploring different subsets of the initial data: involving or not involving the conditional part.

The initial data set comprised a few tautomeric transformations measured at different solvents with significant difference in the equilibrium constant values. These species could be an additional criterion for GTM models quality estimation: thus, their successful separation is an evidence of the model being able to differentiate these objects in spite of the fact that their structures are the same. As a measure of the quality of tautomeric classes separation we applied a special characteristic, Г-score 184 , which can be computed from the GTM class maps and hence does not need the property values be measured in the experiment. The SVR method is used here for the benchmarking purposes.

The article related to this project is given in in the end of this section, with the authorization of all authors.

Modeled object and property

The modeled item is a prototropic tautomeric transformation, for which the equilibrium constant (log K T ) referred to a certain solvent/mixture of solvents and temperature is given. The transformations are assigned to 10 distinct tautomeric classes (Table 5) the active atoms of which, i.e. accepting/donating the hydrogen atom, are attributed and marked. The Marked Atom (MA)-based descriptors are used to compare their performance with the CGR-based approach, that had been tried previously in our group 185 . The structural characterization of the tautomeric transformation was undertaken by describing the "reference" tautomer of each pair. This "reference" tautomer is chosen as the left side of equilibria listed in Table 5, and then coherently applied throughout the work: the "keto" form will be reference for all ketoenol processes. An example of active atoms labeling is given of Figure 21. That is done, at first, in order to avoid repeating description of the same atoms of the second form so that to reduce the overall number of descriptors and, at second, to reduce the number of classes for the depiction and analysis for the case of GTM classification modeling. Thus, the labeling of both active atoms and the preorganized assignation defined by Table 5 allows to fully describe a given equilibrium. The descriptor vector has been composed of the fragments with the labeled donor atom, the fragments with the labeled acceptor atom and the reactions condition descriptors, that were concatenated into a unified descriptor vector representing a single tautomeric transformation of a particular type, measured under a certain reaction condition. Two types of modeling had been performed: the 'structural' and the 'general' ones. The structural models were referred to the GTM classification task solely and were built on a subset of structurally unique transformations, which were 350 out of 695 initial equilibria. The experimental conditions were not included into the descriptor vector for the 'structural' modeling. The 'general' modeling has been performed on the initial data of 695. The descriptor vector in this case included the structures and the conditions. The performance of the GTM for the regression task has been compared with SVR.

Modeling workflow 6.3 Data preparation

The training set, composed of 695 tautomeric transformations, with the values of the logarithm of the equilibrium constant (logK T ) measured in different solvents and at different temperatures, has been critically selected form the database prepared by Gimadiev et al 185 .

The selected dataset contains equilibria for which only two stable tautomeric forms may potentially exist. All transformations are assigned to 10 types of tautomerism (Table 5). The equilibrium constants for them were measures in 12 pure solvents (water, methanol, ethanol, propanol, butanol, cyclohexane, benzene, chlorophorm, DMSO, acetone, DMFA, ethyl ether) and 7 different types of water-organic solvent mixtures (water/ethanol, water/propanol, water/butanol, water/acetone, water/DMFA, water/DMSO, water/ethyl ether) with different proportions of components. The temperatures varied from 233K to 373K.

Type of tautomerism

The number of transformations in the DB

Keto-Enol (I) 271

Amino-Imino (II) 178 Hydrazine-Hydrazone (III) 12

Pyridol-Pyridon (IV) 5

Phenol-Imine -Keto-Amine (V) 33

Thione-Enol -Keto-Thiol (VI) 10

Amine-Thione-Imine-Thiol (VII) 18

Nitro-Aci (VII) 8

Classical Form -Zwitterion (IX) 28

Chain-Ring (X) 132 For some transformations not one, but several different K T values measured at the same conditions were reported in the literature. In this case, logK T for a given equilibrium was calculated as an average of the related experimental values. The structures were standardized by the ChemAxon's Standardizer utility 152 ('basic aromatization' was used).

The subject of unbiased model validation has already been raised in the previous project, where four different scenarios of 'novelty' of the complex with respect to the training set have been introduced (see chapter 5). Regarding to tautomeric transformations, the uncertain/biased statistics could arise for a data set comprising structurally identical equilibria measured in different reaction conditions. In this case, external prediction of the equilibrium constant distinguishes four scenarios: (1) tautomers present in the training set, to be predicted under reaction conditions also seen at training stage (but not in conjunction with those specific tautomers), ( 2 The same equilibria of the initial data set were frequently measured under different reaction conditions, so the number of the unique transformations (i.e. without considering the same structural changes in different solvents as different transformations) in the data base is 350. These transformations have been extracted without the reaction condition part and gathered into the subset of unique transformations, employed in GTM data analysis for the evaluation of the reaction condition influence on data distribution of GTM maps, as well as for the estimation of pure structural clustering as the criteria of appropriateness of a certain descriptor type.

Computational details

The preliminary scanning of 64 different descriptor spaces has been performed by SVR, the MA3 strategy has been used exclusively. The generated structural descriptors have been concatenated with the 14 reaction condition parameters for solvent and temperature (described in section 2.2) and molar fraction of organic component (for water/organic mixtures).

The descriptor set producing the SVR models with the highest R 2 score has been chosen for further evaluation of the remaining labeling strategies. The best descriptor set was based on atom-centered fragments of the length 1-3 and has been used for GTM modeling and visualization and for the external validation of the obtained individual SVR and GTM regression models. The applicability domain method was Fragment Control.

Γ-score. The clustering performance of the GTM can be estimated by the Γ-score 186 which is normalized from 0 to 1 and can be calculated for any data set where the information about classes is available. The Γ-function takes into account k nearest neighbors of each projection.

The more neighbors of each point belong to the same class the higher is the Γ-score. Thus, this score characterizes the quality of class separation on the map. First, for each compound 𝒗 𝒊 , 𝑮(𝒍, 𝒌) should be computed:

𝐺(𝑙, 𝑘) = 1/𝑘 ∑ 𝑔(𝑣 𝑖 , 𝑗) 𝑘 𝑗=1 ( 38
)
where k is the number of nearest neighbors, 𝑔(𝑣 𝑖 , 𝑗) =1 if the jth nearest neighbor of 𝑣 𝑖 belong to the same class, otherwise 𝑔(𝑣 𝑖 , 𝑗) = 0. Then, for each class i, 𝑦 𝑖 (𝑘) is defined as

𝑦 𝑖 (𝑘) = 1/𝑛 𝑖 ∑ 𝐺(𝑙, 𝑘) 𝑛 𝑖 𝑙=1 (39)
where 𝑛 𝑖 is the number of compounds of the class i. And, finally the Γ-score is

Г(𝑘) = 1/𝑁 ∑ 𝑦 𝑖 (𝑘) 𝑁 𝑖=1 ( 40
)
where N is the number of classes. The k number is set to 7.

The details of the computational procedure concerning the specification of the scanned descriptor spaces, the corresponding SVR and GTM parameters of the individual models, as well as the details of data curation and treatment are described in the corresponding sections of the article.

Results and discussions

Data visualization and analysis with GTM

Unique structure subset

The unique data set helps to estimate the influence of the reaction conditions and demonstrates the ability of the GTM method to classify different tautomeric types with and without specification of reaction conditions. The task of predicting the type of tautomeric transformation itself does not have any intrinsic value, because the type of tautomeric transformation can be easily extracted from the transformation equation with a well-known atom-atom mapping without the need to build any models using machine learning methods. Nevertheless, the ability of a given descriptor set to discriminate effectively different types of chemical objects indicates its quality and the ability to be used in building and analysis of different models. The performance of different marked atom strategies has been analyzed with respect to the same descriptor set (the one used herein and after is based on atom-centered fragments, see 6.4).

Figure 23 depicts the GTM classification landscapes for four marked atom strategies for the unique structures subset. Different colors in each landscape correspond to 10 types of tautomerism (Table 5). The corresponding Balanced Accuracy is close to 1 for all tautomeric classes, which correspond to their good separation on the map (App., part III, Table III.3). However, a visual comparison of the landscapes reveals that the labeling strategies MA1 and MA2 separate different classes equally well and better than MA0 and MA3. Thus, the ketoenol and amino-imino classes of tautomeric transformations are well separated from each other and from the other tautomeric types for strategies MA1 and MA2, but not for MA0 and MA3 (Table 6). The corresponding Г-scores for MA1 and MA2 are also higher (Table 6) resulting in maps with more distinct separation of classes compared to MA0 and MA3 maps, where there are more areas of adjoining occupation for different classes.

The entire data set

The entire data set, in which different conditions for the same unique transformation are included in the description, has been visualized and analyzed in the same way as for the unique structure subset above. The corresponding GTM property landscapes, characterizing the distribution of log K T , and the GTM maps colored by classes (class landscapes) for four marked atom strategies are represented in Figure 24.

In comparison with the maps for the subsets of unique structures presented in Figure 23, the maps obtained for the entire data set contain substantially more points, since the individual points on them correspond to different combinations of structural changes with reaction conditions, i.e. different solvents and temperature. The Balanced Accuracy for all tautomeric classes are close to 1 (App., part III, Table III.4), in correspondence with its good separation on the maps, especially for the MA1 and MA2 strategies, the Г-scores of which are higher as well. These results can be explained by the fact that the MA1 and MA2 fragments are more specific for a given type of transformations and therefore different types of equilibria form better separated clusters. The comparison with the unique subset maps shows that the logKT inclusion of the reaction conditions leads to better delimitation of the tautomeric classes: thus, the entire set maps possess less areas occupying by several classes where they are closely located. This could be derived from different solvent affinity and solubility for different chemical transformation, varying structurally within the same tautomeric group, but much more on going from one tautomeric class to another. The GTM method thus clusters data on account of both structures of transformations and the reaction conditions.

Cross-validation of the SVR and GTM models

The performances of the regression SVR and GTM models built on the entire data set for four marked atom strategies are shown in Table 7. According to the results, the MA2 and MA1 strategies lead to the higher predictive performance of GTM models, so as in the case of the SVR models. It could be suggested, that the MA0 strategy is worse because of lack of local descriptors, explicitly accentuating the active atoms and hence bearing the information about the tautomeric type. The weaker performance of MA3 compare to MA1 and MA2 could be due to large number of global descriptors, occulting the influence of both important local descriptors and descriptors of reaction conditions. Notice also that for the GTM models the difference in performances between different strategies is considerably more pronounced than for the SVR models. That could also be due to the fact that the SVR method is known to perform implicit weighting of descriptors, so unimportant global descriptors get low weights, and their adverse effect is minimized. One can also pay attention to the fact that the MA2 strategy very slightly but consistently outperforms MA1 in both GTM and SVR modeling. A putative explanation for this is that the MA2 strategy provides a more detailed description of the environment of the active center due to additional descriptors that are important for modeling of tautomeric transformations. Comparison of R 2 and RMSE values, which characterize the predictive performance of the models, and the Г-scores, which characterize the quality of class (tautomeric type) separation on GTM maps, reveals a consistent correspondence between them: labeling strategies with higher Г-score lead to regression models with higher predictive power. A putative explanation for this is that the map with higher Г-score is characterized by more uniform distribution of data points, which leads to smoother property landscapes with higher predictive power. This opens up interesting prospects for using GTM maps for improving the regression models, because the construction of GTM maps and the maximization of the Г-score for them does not require the knowledge of property values measured in experiment and therefore can be performed for virtual datasets of any size.

GTM solvent separation analysis

The way in which reaction conditions are shaping the logK T property landscapes can be observed by looking at tautomeric equilibria that were studied in different solvents. One can expect that if combinations of the same tautomeric transformation with different solvents are mapped to nearly the same location in the latent space (2D GTM map), while the difference between the values of equilibrium constants for them is significant, then this produces an "activity cliff" which hampers the predictive performance of a models. Otherwise, if they are mapped onto a broad region in the latent space, then the property landscape is smoother, which is favorable for the high predictive performance of the models. To examine this issue, of species on the map, a subset of tautomeric transformations that have a considerable difference (more than 1 in log scale) between the log K T values measured in different solvents have been retrieved. The results are presented in Figure 25 and Table 8. The color of the node corresponds to the property value (log K T ) based on the responsibility contribution of each molecule into the node. According to the obtained picture, the data points are well distinguished on all four GTM maps. The GTM model thus correctly recognizes the dependency of the property from both chemical structure and the solvent nature.

External validation of the SVR and GTM models

To estimate the predictive performance of GTM and SVR models on external test sets, we have chosen the models based on MA2 descriptor set as the one providing the best results. The results of the prediction for the first external set (test set 1), assessing the models performance as a function of experimental conditions, and the second test set (test set №2) assessing the predictive performance for new chemical entities, are given in Figure 26. For both sets the predictive quality of the GTM is better than the one obtained with the SVR model. For the set 2 comprising new structures, thirteen out of 26 transformations were 8) with considerable difference in logKT values due to measurements in different solvents.

found to be out of the model's applicability domain. Statistical parameters calculated for the remaining equilibria within AD show an equivalent performance of both methods (R2 = 0.83 and 0.82; RMSE= 1.0 and 1.2 log units for GTM and SVR, respectively).

Conclusion

This project is devoted to the modeling of tautomeric equilibria accounting for different experimental conditions. At this level, the constituents that need to be considered are the structures of chemical participants, the active centers responsible for chemical transformation and the conditions, under which the latter is conducted. The structures have been encoded with Marked Atom (MA)-based fragment descriptors, where the labels were assigned to atoms donating or accepting hydrogen while the conversion.

Four marked atom strategies have been tried and compared in the frame of the same descriptor type. It was revealed, that the use of labeled (local) fragment descriptors (strategies MA1 and MA2) leads to models with better predictive performance for both, SVR and GTM machine learning methods, rather than the use of unlabeled (global) fragment descriptors (strategy MA0). Moreover, mixing unlabeled and labeled fragment descriptors (strategy MA3) results in a deterioration of the predictive performance in comparison with the use of only the labeled ones (MA1 and MA2). This can be explained by the important role of local structural factors and inessential role of global ones for predicting the constant of tautomeric equilibrium.

Figure 26. Predictive performance of the models on the external sets of new conditions (test set 1) and new structural transformations (test set 2).

The GTM method for the first time has been applied for modeling and visualization of data, which is more complex than the subsets of single molecules ordinarily used for GTM before. The visual GTM analysis, performed using class and property landscapes, has shown that the ability of a descriptor set to provide good separation between different classes of tautomeric transformations correlates with the overall model's predictive performance. For both methods, the best performance has been obtained with MA1 and MA2 strategies, that provides an optimal structural description at the same time not producing a redundant amount of descriptors, that could occult the influence of important local or condition descriptors. The difference in performances for different strategies is less pronounced for SVR models due to the internal mechanism of assigning lower weights to unimportant descriptors in the SVR algorithm.

To quantify the quality of class separation, we applied a special characteristic, Г-score, which can be computed from GTM class landscapes. It was shown, that the descriptor sets providing higher Г-score values in GTM classification landscapes lead to regression models with higher predictive power. This opens up interesting prospects for using GTM maps for improving the predictive models for chemical reactions, since building GTMs having the maximization of Гscore as objective function is a simulation which only requires the knowledge of the reaction types and not of their explicit kinetic or thermodynamic parameters.

The predictive performance of the GTM and SVR models have been compared using two external test sets, providing an unbiased assessment of model's ability in predicting different structures or different reaction condition. The GTM models have shown better predictions (R 2 = 0.62-0.65, RMSE = 0.6-1.96) in comparison with the SVR models (R 2 = 0.39-0.53, RMSE = 0.76-2.2). The GTM approach therefore can be recommended for building other QSPR models, as it combines good predictive performance with the ability to conduct indepth visual analysis of data constituent and of the influence of various factors on quantitative characteristics of chemical processes.

The SVR individual model that includes the assessment of the applicability domain and an automatic labeling of the active atoms, is freely available on our web-server (http://cimm.kpfu.ru/models).

Chapter 7

QSPR modeling and visualization of kinetics properties of cycloaddition reactions.

Cycloaddition (CA) is a classic reaction in organic chemistry being of a fundamental importance for organic synthesis as the main tool for the production of the compounds of cyclic architecture. A large variety of cycloaddition reactions arises from the diversity of reagents, allows to design cyclic adducts of different size, nature and functions. A high regioand stereo-selectivity of cycloaddition explains its wide application at different stages of complex organic synthesis.

In early works, frontier molecular orbitals (FMO) [187][188][189][190][191] calculated by quantum mechanics methods for small congeneric series were widely used for interpretation of thermodynamic and kinetic properties of cycloaddition. Most of modeling studies of the kinetics were conducted using a time-consuming quantum-chemical calculation. An effort has been made to build a linear correlation with experimentally measured physicochemical parameters of the reagents [192][193] . This, however, significantly limits application of the latter to already studied molecules.

Earlier QSPR modeling of chemical reactions were performed on homogeneous series either keeping the solvent, or the structure of the reactants constant [194][195][196][197] . In these studies topological indexes 198 , quantum-chemical 194-195, 199 200 or mixed [201][202][203] descriptors for reagents were used. These works however were restricted to water or gas media exclusively. Such models show high correlation coefficients, but cannot be thought of as universal. For a nonexhaustive overview of the studies carried out in the field of chemical reactivity one can refer to work of A.Warr 204 , I. Baskin et al.,205 or earlier works [196][197]206 .

In spite of significant progress in the field, few attempts have yet been made towards more extensive consideration of chemical reaction: the first one refer to the classification problem in terms of categorical prediction of the conditions (solvent type, catalyst) of the Michael addition 207 . The multicomponent approach of quantitative estimation of the kinetics of chemical reactions explicitly considering both reactants and solvents have been applied to the modeling of the rate constant of S N 2 126, 208 and E2 209-210 reactions leading to a fairly good correlations with the experimental values (R 2 =0.67-0.79). This chapter is devoted to the modeling of kinetic properties of the reactions of cycloadditions comprising (4+2), (3+2) and (2+2) types. The considered properties included the rate constant (logk), the activation energy (Ea) and the pre-exponential factor (logA). Two machine learning methods, SVR and GTM have been involved, where for the latter a unified model, eligible for the prediction, visualization and analysis of all three properties, has been constructed.

Chemical structures of reactants and products were encoded by a single Condensed Graph of Reaction (CGR) representing a given reaction as a single pseudomolecule (described in detail in section 2.1.4). The CGR-based descriptors explicitly characterize the reaction center together with its immediate structural neighborhood. The CGR-based descriptors and the reaction conditions descriptors have been concatenated resulting in descriptors vector describing the entire transformation.

The prepared models have been exhaustively validated on the external test set comprising reactions both structurally different from those in the training set and the same transformations measured under different reaction conditions. Different scenarios of logk assessment were exploited: direct modeling, application of the Arrhenuis equation and temperature-scaled GTM landscapes.

Computational procedure

The modeling workflow is shown in Figure 27. The database of 1849 reactions of cycloadditions, for which the logk, Ea and logA values have been measured in different solvents and at different temperatures, has been collected from the literature (see 7.1.1). The data set has been divided into the training and test set, where the latter has been composed out of 200 reactions picked up randomly from the data set. The remaining 1649 reactions constituting the training set has been used to build three individual models, correspondingly, for logk, Ea and logA prediction. The obtained best SVR and GTM models have been further challenged for an exhaustive external prediction of logk on the prepared test set of 200 reactions.

Data preparation

An initial data set of 2551 reactions of cycloaddition for which all of the reactions contained the experimental measurements of the rate constant (logk), 1356 reactions had activation energies (Ea, in kJ/mol) and 1237 had pre-exponential factor (logA) values was collected manually from the manuscripts of PhD thesis works of Prof. Konovalov's group from Kazan Federal University published in 1970-1990. The data set contained about 85% of Diels-Alder 7.2.3). For the most, the transformations in the test set were measured once, at one temperature. However, there have been 26 reactions with experimental measurements at two different temperatures.

The training set hence encompasses the remain 1649 reactions. The histograms of the distribution for the three properties for the training set is given on Figure 28.

Descriptors

The ISIDA fragment descriptors have been calculated by the in-house ISIDA Fragmentor software [START_REF] Varnek | ISIDA -Platform for virtual screening based on fragment and pharmacophoric descriptors[END_REF] , for which the corresponding CGRs were created from the reaction RDF files using the in-house CGR Designer tool and were stored in modified SDF format. The length of fragments varied for 2 to 14 for sequences and from 2 to 6 for atom-centered fragments. The following options were also used at choice: charges on atoms (Formal Charge), accounting for the terminal atoms of a fragment exclusively (Atom Pairs), and exploring all possible paths instead of shortest paths (DoAllWays). An important option regulating the amount of the overall generated CGR fragments is the 'dynamic bond'. Toggled on, the option produces the fragments, that contains the bonds forming/breaking while chemical reaction and omits the 'generic' fragments, not assigned to the reaction center (see section 2.1.4). That could be used to generate fragments that describe local environment of the reaction center exclusively. For this project, the CGR fragmentation implied the generation of all possible, local and nonlocal, fragments. Overall, 728 descriptor sets have been generated for the preliminary SVR scanning. Structural descriptors have been concatenated with the 14 descriptors of the reaction conditions characterizing solvent and temperature (described in section 2.2).

Building and validation of the models

SVR modeling

SVR models were built and validated using the SVR algorithm implemented in the libSVR package 211 . The modeling was performed using the evolutionary SVR optimizer, 212 which can be used to perform both descriptor space selection and optimization of the operational parameters (epsilon, kernel type, cost, gamma) of the SVR method. The procedure, run for the training set with logk as a modeled property, generated a total of 6000 models. The descriptor set producing SVR individual model of maximal robustness (estimated by the R 2 value) was based on atom-centered fragments of the length 2-3. This descriptor space encompassing 688 descriptors (674 of fragments, 14 of the reaction conditions) was used further for SVR and GTM models building and their external validation.

GTM modeling

The GTM models were built using the evolutionary optimizer 212 which supports the choice of the operational GTM parameters (the number of RBF kernels, the number of grid points, the width factor of radial basis functions, and the regularization coefficient). The descriptor space, based on atom-centered fragments of the length 2-3, chosen by the preliminary SVR scanning was considered. The genetic optimization procedure included 3000 generations. The operational parameters of the GTM method were optimized to predict all three properties, logk, Ea and logA with the highest average R 2 value. The prepared GTM manifold has been used for the visualization of property distribution (property landscapes) and for analysis of separation of three different cycloaddition types on the map (class landscape). The description and the technique of 'coloring' of a GTM manifold with regard to class/property is given in section 3.2.3.

Validation of the models

The performance of the models has been compared by R 2 and RMSE values in 5-fold crossvalidation procedure repeated 10 times after data reshuffling.

Different scenarios of logk assessment for the test set reactions

The logk values for the external set of 200 reactions were obtained using three different approaches: direct assessment, Arrhenius-based assessment and temperature-scaled GTM landscapes. The direct assessment envisages the application of SVR or GTM logk predictive models to the test set reactions.

In the Arrhenius-based assessment the logk values are calculated using the Arrhenius equation (eq. 41). The latter implies the usage of predicted Ea and logA values obtained with the help of the related SVR or GTM models.

ln 𝑘 = ln 𝐴 -𝐸𝑎 𝑅𝑇 ⁄ (41) Rate constant is more sensitive to temperature compared to any other individual structural or solvent descriptors. This effect could hardly be accounted for by GTM-based model where each descriptor has exactly the same weight. Therefore, the temperature-scaled approach has been proposed for this purpose. It implies the construction of series of logk GTM landscapes each corresponding to a specific temperature range.

Scaling of logk measured at temperature T1 to temperature T2 can be performed according to Van't Hoff equation: 𝑙𝑜𝑔 The entire temperature range (273K-423K) has been divided into 15 subranges of 10K each, e.g. from 273K to 283K. For each subrange, the corresponding logk landscape was constructed, i.e., each logk was recalculated to median temperature of a subrange (e.g., 278K for the range 273K-283K) using Van't Hoff equation. It should be noted, that the error related to temperature deviation of 5K (0.15-0.3 logk unit) is far less than the related RMSE values (Table 10). Each of 15 temperature-scaled logk landscapes were calculated from experimental logk values of the training set using the Van't Hoff relationship, log 𝑦 = 0.26 and related median subrange temperature. At the validation stage, the program selects particular logk landscape that will be used for prediction, corresponding to the external reaction temperature.

Results and discussions

GTM visualization of the training set

Generative Topographic Map shown on Figure 29 shows well separated zones populated by the (4+2), (3+2) and (2+2) cycloaddition reactions. The ratio of sizes of these zones correspond to classes occurrences in the training set: thus, the major class (4+2) occupies the largest area, the following after is the (3+2) zone, which in turn is slightly bigger than the zone populated by (2+2) reactions. A distinct separation of cycloaddition classes on the map infers a competency of the chosen descriptors space, able to correctly discriminate reactions belonging to different types.

Another technique of GTM visualization is the property landscapes, characterizing the training set distribution of a certain property. Figure 30 shows GTM property landscapes for logk, Ea and logA, where the gradient from blue to red corresponds to the property value variating from low to high (in the range its own for each property). Simultaneous analysis of all three landscapes helps to understand decomposition of logk on Arrhenius equation parameters logA and Ea (eq. 41). According to collision theory, A is the frequency of collisions in the "correct orientation". Thus, logA values might be related to steric interaction of reactants in a prereaction complex. Activation barrier Ea accounts for electronic and steric effects in the transition state. Typically, high reaction rates correspond to large logA and small Ea (zone A on the maps) whereas low logk values result from a combination of small logA and large Ea (zone B). However, in some cases low activation barriers doesn't compensate too low logA (zone C), or, on the contrary, large logA is associated with high activation barriers (zone D).

Zones Example of reaction logk

Cross validation of the SVR and GTM models

The cross-validation result (5CV) for all three properties, logk, Ea (kJ/mol) and logA, for SVR and GTM models are shown in Table 10. The SVR modeling involved the construction of three different models for each property, whereas for the GTM model predictions were assessed with three property landscapes built on one same manifold (see section 7.1.3). The performances of SVR models exceed those of GTM model. However, the overall statistics for GTM looks reasonable for all three properties.

Property

External validation of the SVR and GTM models

The external test set has been predicted in three different ways (see 7.1.3), two common for both, SVR and GTM are (i) the direct logk prediction and (ii) the prediction of the logk using the Arrhenius equation (eq. 41) and the predicted values of Ea and logA. The results for both strategies are given in Table 11 and Figure 31. The weaker results for the Arrhenius-derived logk prediction could be explained, first, by the accumulation of errors and, second, by the narrower Applicability Domain (AD). The accumulation error comes from the prediction of a property by means of the accessory predicted values. The many are the latter, that bigger is the chance of misprediction, affecting the overall accuracy. The second reason is related to the fact, that the amount of experimental data for Ea and logA is smaller by one third than for logk, meaning less accurate training and a narrower AD. The corresponding statistics accounting for AD (Fragment Control) for the Arrhenius-based prediction is much better, where the number of the compounds is however shrinked from 200 to 164. Temperature-scaled logk GTM landscape. Reaction rate assessment with temperature-scaled logk landscapes shows a minor predominance over the direct GTM logk prediction (Table 11). This trend was confirmed in logk predictions for a set of 13 reactions, each measured at two different temperatures (Figure 32). The temperature-scaled approach is less prone to return the same rate constant values at different temperatures, showing an overall better reproduction of temperature dependence. Direct log k assessment

Conclusion

This project is devoted to the modeling of reactions of cycloaddition, with structurally varying reagents, measured at different reaction conditions. Involvement of multiple atoms and bonds during the chemical transformation needs an efficient way of structural description. Here the structures have been encoded with the Condensed Graph of Reaction (CGR), that explicitly designates the bonds forming/breaking while the reaction transformation. The data set of 1849 reactions, comprised of (4+2), (3+2) and (2+2) types of cycloaddition, has been attributed with three experimentally measured properties: rate constants (logk), activation energies (Ea, kJ/mol) pre-exponential coefficients (logA).

Demonstrated a strong ability in the managing of structural/condition chemical data in the previous project of tautomeric equilibria modeling, the GTM method is challenged here in prediction of several kinetic properties with one single GTM model. The SVR models have been built individually for each property. The cross-validation R 2 range for the GTM model is 0.62-0.80 and for the SVR models is 0.89-0.94.

The external validation has been performed on a set of 200 molecules, randomly picked from the initial dataset. Two common ways of logk estimation have been fulfilled: the direct prediction with the corresponding logk models, and the Arrhenius-based recalculation through the predicted Ea and logA values. For both methods, the direct prediction is more accurate than the Arrhenius-based assessment (RMSE direct = 0.5 (SVR) and 0.9 (GTM); RMSE Arrhenius = 0.93 (SVR) and 1.23 (GTM)), which is mostly related to the fact that the amount of data for logA and Ea is much less, than for the logk training set, meaning narrower AD. Thus, Fragment Control AD significantly improved the performance of both SVR and GTM models.

The overall results of SVR method, having the mechanism of implicit ranging of descriptors according to their impact, are better, than for a single unified GTM model, where all the descriptors are taken into relatively equal account. To enhance the predictive ability of GTM with respect to temperature, the temperature-scaled landscapes have been constructed. The performance of the approach was the same as for the direct logk GTM assessment, with a slight improvement in temperature dependency reproduction with regard to 13 reactions measured at two different temperatures (RMSE direct = 0.48, RMSE temperature-scaled = 0.39).

Using a single GTM map provides another advantage: thus, since the distribution of the reactions are the same for each property, it is possible to analyze the dependencies of the properties and the structural features they could be determined by. Thus, simultaneous analysis of property landscapes for logk, logA and Ea helped to identify different types of reactions with respect to the interplay between the Arrhenius parameters logA and Ea. With respect to class separation, the three classes, (4+2), (3+2) and (2+2) are well-separated on the map, occupying the areas proportional to their share in the training set.

The SVR individual models for the rate constant, activation energy and pre-exponential coefficient, supporting the choice of solvent and temperature, are freely available on our webserver: http://cimm.kpfu.ru/models.

Chapter 8

QSPR modeling of the rate constant of S N 1 reactions.

In this chapter we report predictive models for the rate constant of S N 1 reactions involving both Marked Atom-based and the Condensed Graph of Reaction-based approaches. Unimolecular nucleophilic substitution (S N 1) is a reaction with a two-stage mechanism, in which the heterolytic bond cleavage of a neutral molecule precedes the reaction with the nucleophile (Figure 37). Since the rate-determining step of a bond cleavage involves only a substrate, S N 1 is the first-order reaction.

The factors affecting the rate constant of S N 1 reactions were investigated in numerous quantum mechanical calculations of different levels (HF, MP2, G2) [211][212][213][214][215] . These works are useful in investigation of various factors, such as electronic effects of substituent or size of the molecule, influencing the rate constant. The molecular dynamics allows to go further and investigate the reaction process altogether with solvent effects: these works 216-217 help to evaluate the effects of the solvation energy (water only) to the feasibility of the reaction process. Some of the experimental studies were devoted to the substituent effects 192 or certain steric influence 218 (ortho-effect). However, all these studies were performed on small series of reactions and did not account for solvent effects, or only water medium was considered.

In order to generalize the substituent effects on the rate constant, several scales have been proposed. Thus, the α constants introduced by Uggerud 219 has been designed to rationalize the substituent effects of the alkyl groups. A comprehensive nucleofugality scale has been proposed by Mayr [220][221]222 . Later on, the applicability of this scale was expanded for the prognosis of the preferable, S N 1 or Sn2, mechanism 223 . It should be noted, however, that the borderline between S N 1 and S N 2 mechanisms is the subject of considerable controversy. Thus, opposed to Ingold 224 , considering the two mechanism to be distinct, discrete processes it has been ascertained that the clear-cut distinction for many cases could not be established due to gradual transformation of one mechanism into the other. The borderline cases of the concurrent S N 1 -S N 2 reactions have been pointed out for benzylation of pyridines 214 , benzophenones and benzhydrols 223 , arylbromoethanes 225 , tosylates 226-227 and 4-methoxybenzyl derivatives 228 . To interpret these interjacent cases, Winstein's and Sneen developed a concept of different types of ion pairs intermediates [229][230][231][232] . Schleyer and Bentley criticized this concept and suggested that there is a gradation of transition states between the S N 1 and S N 2 extremes [233][234] . The differentiation of the mechanism, however, to a fair degree of accuracy could be delineated by the lifetimes of the potential intermediates [235][236] . The abovementioned scales could be a useful tool for reactivity analysis only within a congeneric series of substrates reacting under similar conditions. The first attempt to develop a more generalized model was performed by Kravtsov et al. 237 on a set of 1661 reactions studied in different solvents. The model involving fragmental descriptors and Fukui indices for structures and solvent descriptors performed well in cross-validation (R 2 = 0.75 and RMSE = 0.61 logk units).

The goal of this work is to build a model for the rate constant of S N 1 reactions applicable for a wide variety of structures and accounting for reaction conditions. A large data set of 8056 reactions, representing the first dissociation step of S N 1 reactions (Figure 37), was used in SVR model building. The molecules have been classified according to the atoms attached to the cleaved bond (C-Hal, C-C, etc). The GTM method was employed for the purposes of data visualization. The most robust SVR model was applied to two external test sets, evaluating the model performance to predict new structures and new reaction conditions.

Computational procedure

The workflow of the rate constant modeling is given in Figure 34. The training set of 8056 reactions of dissociation, has been encoded by CGR-based and MA-based fragment descriptors. The best descriptor space for MA and CGR approaches, producing the most robust SVR models has been chosen for the prediction of two external test sets.

Data preparation

The data set of 11748 transformation measured in 28 different solvents under the temperature range from 202K to 528K, has been collected form the literature [START_REF] Palm | Fundamentals of Quantitative Theory of Organic Reactions[END_REF]170 . For each of the reaction, the experimental rate constant value (logk) has been attributed. The transformations have been represented as the first step of S N 1 reaction, the heterolytic bond dissociation (Figure 37). . The data set underwent an exhaustive cleaning protocol of removing duplicates, inorganic and wrong-drawn compounds. For more nearly 2000 compounds, the reaction rate has been assigned to the catalytic reaction and for some of the cases the catalyst (carbonic acids' salts) significantly influenced the reaction rate.

As the structure of the catalyst is not included into consideration, these reactions have been removed. The data set contained a lot of duplicates, when the reactions are the same in terms of both, structure and condition, but the experimental rate constant for them were different. The amount of these species with respect to the corresponding difference in log k values is shown in Figure 35. The molecules of the first series (the logk difference is within 0.1-1) that had same structures and same reaction conditions have been merged and an average logk value has been assigned. For the cases with logk difference more than 1, all the duplicates have been removed.

The prepared data set of 8456 heterolytic dissociations could be classified in five classes, upon the basis of the type of breaking bond,correspondingly:

•

Descriptors

Both types of ISIDA descriptors, MA-based and the CGR-based, have been used (described in section 2. 1.4). For the MA-based descriptors, the atom neighboring to the leaving group was marked only (Figure 37). fragments. The Formal Charge, Atom Pairs and DoAllWays options were also used. The structural descriptors have been concatenated with 14 reaction condition parameters reflecting solvent and temperature (described in section 2.2). Overall, 270 CGR-and MA-based descriptor spaces have been generated and examined.

Building and validation of the models

SVR modeling

SVR models were built and validated using the SVR algorithm implemented in the libSVR package 211 . The modeling was performed using the evolutionary SVR optimizer 212 . The procedure run for the training set and generated a total of 3000 models. The best descriptor spaces for each descriptor type (MA or CGR), producing the SVR individual models of maximal robustness (selected descriptor space, prescribed kernel type, epsilon, gamma and cost parameters), have been chosen for the cross-validation comparison performance and the external test sets prediction. The winning CGR descriptor set is based on sequences of atom and bonds of the length from 1 to 5, accounting for the charge of the atom ('Formal Charge option) and encompasses 1186 descriptors overall. The winning MA descriptor set is based on atom-centered fragments of atoms and bonds of the length 1-4, accounts for atom's charges and encompasses 996 descriptors.

Validation of the models

The performances of the two SVR models have been compared by R 2 and RMSE values in 5fold cross-validation procedure repeated 10 times after data reshuffling and on external validation on two external test sets.

Results and discussions

GTM visualization of the training set

The GTM method has been employed for data visualization. The training set reactions were classified according to the type of cleaved heterolytic bond: C-Hal, C-O, C-C, C-S or S-O (S-N). Since experimental conditions had no impact on this analysis, a subset of 1820 reactions differed only by their structures were considered. The GTM manifold has been built on the CGR-based descriptors with default parameters. Figure 38 shows that all 5 classes of reactions are well separated on the map.

Cross validation of the SVR models

The performance of the two winning SVR models are represented in Table 13.The CGRbased descriptors explicitly describe the reaction center, corresponding to the breaking bond.

The MA-based descriptors took into account an atom neighboring to the leaving group. It means that for the MA-based descriptors, the structure of the leaving group was encoded, but it's association with the active center was not specified. Similar performance of both models suggests that both strategies of reaction center encoding provide with similar description of a reaction.

Descriptors

External validation of the SVR models

The results of the external validation are shown in Table 14 and Figure 39. The predictions for reactions already occurred in the training set (test set 1) are slightly better than for structurally new reactions (test set 2). The performances of both descriptor types are similar, as it was observed on the cross-validation stage.

That should also be noted, that for some of the reactions a concurrent S N 1-S N 2 mechanism could emerge (examples are given in 'Introduction' of this section), thereby the experimental rate constant would be affected by the competitive mechanism, that would lead to experimental uncertainties. Furthermore, the S N 1 mechanism is often accompanied with the rival E1 reaction, that could complicate the retrieving of 'pure' S N 1 rate constants. The mentioned facts should be taken into account when reckoning the expected predictive accuracy of any S N 1 model.

Table 14. Predictive performance of the SVR models on the external test sets of new conditions (test set №1) and new structural transformations (test set №2

). The predicted property is the Sn1 rate constant.

Conclusion

A large dataset of 8056 S N 1 reactions proceeding in various solvents and at different temperature has been collected from the literature. The data was visualized with the help of Generative Topographic Mapping. Reactions with different types of reaction centers were well separated on the map.

Predictive SVR models for the rate constants of S N 1 reactions were built using two types of descriptors: Condensed Graph of Reaction-based and Marked Atoms-based. Both types of models performed well in cross-validation (RMSE = 0.51-0.52) and on the external test sets. The model predicts logk better for the reactions with new conditions (RMSE = 0.68-0.70) than for the reactions with new structures (RMSE = 0.87-0.90).

The model is suitable for the prediction of both mono-and polyfunctional species. For the latter, the model should be provided with an input file (.sdf) with labeled active centers, one label per one molecule. Thus, a molecule with two binding centers should be written in the input file twice, with the first labeled active atom and with the second.

Tautomeric equilibria of different tautomeric classes.

The model predicts tautomeric equilibrium constant (log K T ) and is available on the webserver https://cimm.kpfu.ru/. The SVR consensus model consists of ten best individual models with R 2 ranging from 0.75 to 0.82. The employed descriptors spaces are listed in Appendix (part III, Table III.5). Ten tautomeric classes are supported (see Table 5). The prediction is rendered for the left-to-right type of equilibria, as listed in Table 5, and not vice versa. That means that the user is expected to write ketone transforming to enol to obtain a keto-enol equilibrium constant, the enol-ketone query would not pass. The models are based on Marked Atom descriptors. The labeling of the corresponding atoms is done automatically. The model supports:

• Various reaction conditions: solvent, mixture of solvents, temperature (Figure 41, a)

Kinetic properties of cycloaddition reactions

Three SVR individual models devoted to the kinetic properties of reactions of (4+2), (3+2) and (2+2) cycloadditions are available on the web server https://cimm.kpfu.ru/. The models are built on atom-centered CGR-based fragments of the length 2-3. Three properties are considered: the rate constant (logk), the activation energy (Ea) and the pre-exponential factor (logA). The performances of the models are given in Table 10. For each of the model, variability of the reaction conditions (solvent, temperature) (Figure 42, a) as well as the estimation of the level of trust of the prediction are supported (Figure 42,b).

Conclusion

This thesis has been devoted to the modeling and visualization of chemical interactions with the aid of local descriptors, identifying and explicitly distinguishing the dynamic nature of molecular sites, directly involved into a chemical process. The consideration has been carried out with respect to the sophistication of the modeled chemical object, correspondingly from intermolecular interactions to chemical reactions, varying reagents and different reaction conditions of which have been explicitly taken into account.

The developed methodology represents a chemical interaction by local fragment descriptors, encoding structural features of the interacting molecules, coupled with special physicochemical parameters of the experimental conditions (solvent, solvent mixtures, temperature). The Marked Atom (MA)-based and the Condensed Graph of Reaction (CGR)based ISIDA fragment descriptors have been employed. The MA-based local descriptors have been chosen for the description of chemical processes, the active site of which involved no more than two atoms, since otherwise, the length of the descriptor vector could be too large. These were the cases of intermolecular interactions and tautomeric equilibria. For the case of chemical reactions, involving different reactants and undergoing multiple bond cleavageformation, the CGR-based fragments, enable to efficiently encode chemical transformation in a condensed, concise form, have been used. The local approach has been successfully used to predict different thermodynamic and kinetic properties of halogen and hydrogen bonding, tautomeric equilibria, cycloaddition and S N 1 reactions. The accuracy of the developed QSPR models are close to experimental errors.

It has been demonstrated, that the models, trained on the complexes with single halogen or hydrogen bonds were able to predict stabilities of the complexes with multiple bonds of these types. This opens a perspective to use the models in computer-aided drug-and supramolecular systems design.

For the first time, the GTM method has been employed to model and visualize entire chemical processes instead of individual species. Thus, on the example of tautomeric equilibria it was shown that the species, measured in different solvents are well separated on the map, meaning that the GTM-based model could successfully perceive the dependence of the modeled property as a function of two equally important components, the structural and the reaction conditions one. The capability of GTM in building one single model able to predict different kinetic properties of chemical reactions has been demonstrated on the example of cycloaddition. For the classification task, employed in the projects devoted to tautomeric equilibria, cycloaddition and S N 1 reactions, a good class-determination performance, resulting in a perfect visual separation of the classes on the GTM map, has been accomplished.

The developed QSPR models for the prediction of the stability of halogen-bonded complexes with I 2 , of the equilibrium constant of tautomeric transformations and of the rate constant, activation energy and the pre-exponential factor of cycloaddition reactions are publicly available via our internet-platforms https://cimm.kpfu.ru/ and http://infochim.ustrasbg.fr/webserv/VSEngine.html. The implementation incorporates an automatic detection and labeling of the active atoms, or CGR generation, as well as the applicability domain estimation. -6.45 -7.38 -8.85 -14.83 -11.09 -13.30 -11.58 -10.89 -13.07

 G pred = 0.6* Σ G pred,i IIRB-MA3--FC-1-4 -s 3 -t 1 -c 0.01005184 -e 1.208808 -g 1.041811 -r 6.9

Descriptor set LibSVM setup 5CV R 2 5CV RMSE

IAB-MA2--P-2-14 -s 3 -t 2 -c 2.446919e+02 -p 1.576160e-01 -g

Tableau 1 .

 1 hydrogène. Dans ce second exemple, les modèles sont construits sur des données obtenues en environnement contrôlé : un solvant de référence (CCl 4 ) et une température (298K). Mais cette propriété fait maintenant intervenir deux espèces chimiques : l'accepteur et le donneur de liaison hydrogène. Ceux-ci doivent être simultanément pris en compte dans les modèles et dans l'évaluation des performances de ces modèles. Le troisième exemple illustre un niveau de complexité encore supérieur : la modélisation des équilibres tautomères. Une molécule peut exister sous plusieurs formes qui ne se distinguent que par la position dans la structure chimique, d'un ou plusieurs atomes d'hydrogènes. Ces différents états sont les tautomères et leur prévalence relative est contrôlée par une constante d'équilibre tautomères. Cette constante dépend, en fait, de conditions expérimentales : solvant et température. Afin d'en tenir compte explicitement, il est nécessaire d'introduire des variables supplémentaires et pertinentes pour les décrire. Le dernier exemple traite des réactions de cycloaddition et de S N 1, où différents réactifs et différentes conditions réactionnelles sont impliqués (Tableau 1). Ce projet met en oeuvre tous les développements présentés auparavant. Informations sur les modèles prédictifs développés dans ce travail. sur des atomes marqués (MA) b descripteurs basés sur graphes condensés de réaction (CGR)

1 . 1

 11 thèse comprend sept chapitres. Le premier chapitre d'introduction décrit divers descripteurs locaux utilisés dans la modélisation. Le deuxième chapitre fournit des informations sur les outils chémoinformatiques utilisés dans cette étude. Les chapitres 3 et 4 décrivent des modèles prédictifs pour évaluer les stabilités des liaisons halogènes et hydrogènes, respectivement. Les chapitres 5, 6 et 7 sont consacrés à la modélisation prédictive de certaines propriétés thermodynamiques et cinétiques des réactions chimiques : la constante d'équilibre tautomère et les constantes de vitesse de la cycloaddition et des réactions S N 1. Descripteurs locaux ISIDA pour la modélisation, l'analyse et la visualisation de données chimiques Les descripteurs locaux utilisés dans ce travail sont des sous-ensembles des descripteurs ISIDA 1 représentant des fragments de différentes longueurs et topologies d'un graphe moléculaire donné. Ces fragments contiennent au moins un atome ou une liaison étiqueté. Deux types de fragments ont été considérés: basés sur les atomes marqués et sur le graphe condensé de réaction (voir la Figure 1). Les atomes marqués (MA) sont des atomes d'une structure chimique, annotés pour leur pertinence vis-à-vis d'un problème donné. Dans ce travail, ce sont les atomes impliqués dans les interactions intermoléculaires. En d'autres termes, ce sont des atomes donneurs d'électrons dans des accepteurs de liaison halogène, donneurs et accepteurs de proton dans les espèces formant des liaisons hydrogènes, les atomes qui perdent / reçoivent des protons dans des équilibres tautomères.

Figure 1 .

 1 Figure 1. Exemple de structures pour lesquelles des descripteurs locaux basés sur MA (en haut) ou basés sur CGR (en bas) ont été générés. En haut: dans le complexe lié à l'hydrogène, les croix désignent les Atomes Marqués. En bas: dans le Graphique Condensé codant pour la réaction de cycloaddition (4 + 2), les points et les tirets représentent respectivement des liaisons chimiques formées et rompues. Quelques exemples de descripteurs générés sont donnés sur la droite.

Figure 2 .

 2 Figure 2. Valeurs de logK BI2 prédites vs expérimentales pour l'ensemble de test externe (à gauche) et des exemples d'évaluation de logK BI2 pour deux molécules polyfonctionnelles (à droite).

Figure 4 .

 4 Figure 4. Les énergies libres prédites vs expérimentales (kJ / mol) pour les complexes 1:1 avec deux liaisons hydrogènes coopératives (à gauche) et les deux exemples de ces complexes (à droite).

  Figure 3. Les énergies libres prédites vs expérimentales (kJ / mol) pour l'ensemble de test externe de stabilité de la liaison hydrogène mesurés dans différents solvants.

Figure 5 .

 5 Figure 5. Valeurs de constante d'équilibres tautomères pour l'ensemble de transformations étudié dans de nouvelles conditions (à gauche) et l'ensemble de test contenant de nouvelles transformations (à droite).

  ) et (2+2), associées à leurs valeurs expérimentales de la constante de vitesse (logk), à 1356 valeurs des énergies d'activation (Ea) et à 1237 valeurs du coefficient préexponentiel (logA). Les descripteurs basés sur les graphes condensés de réactions ont été utilisés pour construire des modèles SVM et GTM individuels, suivis de leur validation (les valeurs du log k ont été prédites) sur le jeu de test de 200 réactions sélectionnées aléatoirement dans la base de données. Les énergies d'activation et le facteur pré-exponentiel ont été modélisés. Ceux-ci sont utilisés pour estimer la constante de vitesse, logk, en suivant une loi d'Arrhenius. La Figure 7 représente les paysages d'activité généré par la GTM, pour le logk, l'Ea et le logA. Les résultats sont en accord avec des concepts chimiques généraux, les réactions à faible logk, projetées dans les zones de faibles logA bas et Ea sont caractérisées par d'importantes contraintes stériques (Figure 7, a). Pendant ce

Figure 6 .Figure 7 .

 67 Figure 6. Paysage de classe GTM (à gauche) et paysage d'activité (au milieu). Les réactions 1 (a, b) et 2 (a, b) sont données à droite. Les nombres entre parenthèses correspondent aux valeurs logK T .

2. 5

 5 Modélisation QSPR de la constante de vitesse des réactions S N 1.Un ensemble de 8056 données de réactions S N 1 a été utilisé pour construire des modèles de régression SVM de leur vitesse de réaction. Les réactions étaient codées par des descripteurs intégrant des atomes marqués de type MA3 ou calculés sur de structures CGR. Les modèles ont été validés sur deux jeux de données supplémentaires: l'un contenant des réactions connues étudiées dans de nouvelles conditions (test 1) et l'autre de nouvelles réactions (test 2). Le modèle fonctionne de manière similaire sur les deux ensembles de test: R 2 test1 = 0.64-0.67, R 2 test2 = 0.55-0.58, RMSE test1 = 0.68-0.70, RMSE test2 = 0.87-0.90.

Figure 8 .

 8 Figure 8. Modélisation des logarithmes de constantes de vitesse : valeurs prédites vs valeurs expérimentales pour le jeu de test (à gauche). Séparation de trois classes de réactions (4+2, 3+2, 2+2) en utilisant la méthode GTM (à droite).

Figure 5 .

 5 Figure 5. Modélisation des logarithmes de constantes de vitesse : valeurs prédites vs valeurs expérimentales pour le test externe étudiées dans de nouvelles conditions (à gauche) et de test externe contenant de nouvelles transformations (à droite).

▪▪

  ) la sélection et l'élaboration d'une combinaison optimale de descripteurs caractérisant la structure chimique, d'une part, et les conditions expérimentales, d'autre part; (3) la capacité de la GTM à visualiser et à modéliser des processus chimiques entiers (structures et conditions); (4) la capacité des modèles formés sur les complexes avec une seule liaison halogène ou hydrogène à prédire la stabilité des complexes avec de multiples liaisons de ces types. Une combinaison de descripteurs de fragments locaux décrivant des structures moléculaires et de descripteurs spéciaux caractérisant les conditions expérimentales (solvant et température) a été utilisée avec succès pour développer des modèles prédictifs de certains paramètres cinétiques et thermodynamiques des liaisons halogènes et hydrogènes, des équilibres tautomères et de deux types de réactions chimiques. Les descripteurs les plus appropriés pour les réactions chimiques combinant plusieurs réactifs et produits sont ceux générés à partir des Graphes Condensés de Réaction. Sinon, diverses stratégies utilisant les Atomes Marqués ont été recommandées. ▪ Les modèles construit sur des mesures réalisées sur des complexes impliquant une seule liaison halogène ou une seule liaison hydrogène ont pu être utiliser avec succès pour estimer les stabilités des complexes contenant plusieurs liaisons de ces types. Cela ouvre des perspectives pour utiliser ces modèles dans la conception assistée par ordinateur de nouveaux systèmes supramoléculaires. Pour la première fois, la méthode GTM a permis de visualiser des processus chimiques décrits par l'ensemble de leurs réactifs et de leurs produits et leurs conditions réactionnelles plutôt que par des espèces chimiques individuelles. Ainsi, sur l'exemple des équilibres tautomères, il a été montré que les espèces mesurées dans différents solvants sont bien séparées sur la carte. Dans l'exemple de la cycloaddition, les trois types de réaction sont bien séparés aussi. ▪ Les modèles QSPR développés prédisant les constantes d'équilibre des transformations tautomères, la stabilité des liaisons halogène avec I 2 et la constante de vitesse, le facteur pré-exponentiel et les énergies d'activation des réactions de cycloadditions sont disponibles pour les utilisateurs via nos plateformes internet https://cimm.kpfu.ru/ et http://infochim.ustrasbg.fr/webserv/VSEngine.html.

Figure 1 .

 1 Figure 1. Transition state for the acid (a) and base (b) hydrolysis of esters.

Figure 2 .

 2 Figure 2. The structure of bicyclo-[2.2.2]-octane-l-carboxylic acids.

Figure 3 .

 3 Figure 3. Example of structures for which MA-based (top) or CGR-based (bottom) local fragment descriptors were generated. Top: in the hydrogen-bonded complex, the stars denote Marked Atoms. Bottom: in the Condensed Graph encoding the (4+2) cycloaddition reaction, dots and dashes represent, respectively, formed and broken chemical bonds. Some examples of generated descriptors are given on the right.

Figure 4 .

 4 Figure 4. Examples of ISIDA MA descriptors (sequences of the length 4) generated for different marked atom strategies.

Figure 6 .Figure 5 .

 65 Figure 6. An example of ISIDA atomic logP increments. The atoms are divided into groups, which are denoted by the color code, in accordance to their calculated logP increments.

  Figure 7. Hierarchical clustering dendrogram based on the 13 solvent parameters used for the solvent effect account in the application part of the thesis.

Figure 8 .

 8 Figure 8. An illustration of the SVR algorithm: a separating hyperplane (denoted with red line) is drawn so that to separate the objects belong to different classes with the maximal distance betweenthe nearest points and the hyperplane.

Figure 9 .

 9 Figure 9. The kernel trick of the SVR: a nonlinear transformation from feature space to higher-dimensional space where the objects are easily separated.

Figure 10 .

 10 Figure 10. An illustration of how the data points are fitted by the GTM manifold (left)and further projected on the GTM map (right).

Figure 11 .

 11 Figure 11. Possibilities of GTM visualization: the class landscape (left), representing the separation of 10 different tautomeric classes, and the property landscape (right) characterizing the distribution of the equilibrium constant values over the map.

Figure 12 .

 12 Figure 12. Schematic representation of 5-fold cross-validation procedure. Initial dataset is divided into 5 parts, on each fold a model is trained on 4 parts and is applied to predict the last one. At final, all predictedvalues are gathered for statistical evaluation.

Figure 14 .

 14 Figure 14. The electrostatic potential surfaces of thrifluoromethyl halides where the areas of positive electrostatic potential correspond the 'σ-holes' that determine the capability and the activity in halogen-bond interaction.

Figure 16 .

 16 Figure 16. Predicted vs experimental log KBeff for the external test set of 11 polyfunctional molecules.

Figure 17 .

 17 Figure 17. Workflow of the modeling of the Gibbs Energy (ΔG, kJ/mol) of hydrogen bonded complexes of different donors and acceptors.

  were given for different solvents, meaning that the predicted values of ΔG underwent the solvent-specific corrections. The corrections have been obtained with the help of linear solvents correlations, retrieved from the initial data, where some of the complexes, apart from the values in CCl 4 , had additional measurements in other solvents. The second external set (test set 2) contains 12 dimers with cooperative HBs (Figure18) measured at 'standard' (meaning the same as for the training set) conditions.

Figure 18 .

 18 Figure 18. An example of complex with two cooperative hydrogen bonds.

Figure 21 .

 21 Figure 21. Example of atom labeling to encode a tautomeric process. The right-hand side tautomer results from the motion of proton from donor atom (red star) to acceptor atom (green star) in the left-hand side tautomer.

Figure 22 .

 22 Figure 22. Workflow of the modeling of tautomeric equilibria constant (log KT).

Figure 23 .

 23 Figure 23. GTM classification landscapes for four marked atom strategies for the unique structures subset of 350 equilibria.

Figure 24 .

 24 Figure 24. GTM property (top) landscapes, representing the distribution of the logKT values, and the GTM classification (bottom) landscapes, representing the separation of 10 different tautomeric classes, built for the entire data set of 697 equilibria. The descriptors spaces for four marked atom strategies are based on atomcentered fragments of the length 1-3.

Figure 25 .

 25 Figure 25. The solvent separation of structurally same pairs of tautomeric equilibria (Table 8) with considerabledifference in logKT values due to measurements in different solvents.

Figure 27 .

 27 Figure 27. Workflow for the modeling of the rate constant (log k), activation energy (Ea) and pre-exponential coefficient (log A) of the reactions of cycloaddition.

Figure 28 .

 28 Figure 28. Distribution of rate constant, activation energy and pre-exponential factor values in the training set of 1649 cycloaddition reactions.

  The temperature coefficient log 𝑦 was computed as an average of log 𝑦 𝑖 (i = 1 -n) values calculated for n series of reactions, each studied at several temperatures. Overall, the reaction rates of 358 reactions in the training set ware measured at 2 -6 different temperatures. The smallest temperature difference in a series was 10K and the largest was 110K. The temperature coefficient log 𝑦 varied from 0.08 to 0.59; its average value was 0.26.

Figure 29 .

 29 Figure 29. GTM class landscape displaying separation of three different types of cycloaddition in the training set of 1649 reactions.

Figure 30 .

 30 Figure 30. GTM property landscapes for the rate constant (logk), activation energy (Ea) and the pre-exponential coefficient (log A) for the training set of 1649 reactions of cycloaddition. The encircled areas correspond to different combinations of logA and Ea contributions into logk.

Figure 31 .Figure 32 .

 3132 Figure 31. Predicted vs experimental log k values obtained with the direct log k prediction (left) or with the Arrhenius-based recalculation (right) for the external test set of 200 reactions of cycloaddition. The statistics is given inTable 11.

Figure 34 .

 34 Figure 34. Workflow of the modeling of the rate constant of the dissociation step of Sn1 reactions.

  class 1: C-Hal (Hal=F,Cl,Br,I), ~50% of the data • class 2: C-O, ~45% of the data • class 3: C-C, ~3% of the data • class 4: C-S, ~1,5% of the data • class 5: S-O and S-N, ~0.5% of the data This classification will be further used for the GTM visualization. The data set has been divided into the training and test set data. Two external sets, each containing 200 dissociations have been prepared. The first one consists of the transformations that have been already occurred in the training set, but under different reaction conditions. This set is aimed at assessing the model's ability to predict the rate constant in new reaction conditions. The second test set comprised the structures which are different from the ones in the training set and thus assessing the model's ability to manage with structurally new transformations. The training set included the remaining 8056 reactions. The histograms of logk distribution of the training set is given in Figure 36.

Figure 35 .

 35 Figure 35. Occurrence frequency distribution histogram indication the amount of duplicates with a certain (log unit) difference in log k values.

  The CGR-based fragments have been counted by the in-house ISIDA Fragmentor software[START_REF] Varnek | ISIDA -Platform for virtual screening based on fragment and pharmacophoric descriptors[END_REF] , for which the corresponding CGRs were created from the reaction RDF files using the inhouse CGR Designer tool and were stored in modified SDF format. The labeling of the structure for MA descriptors was performed through the CGR, the 'dynamic charge' option of which allows unambiguously locate the atom acquiring positive or negative charge during the reaction. The degrees of 'locality' and the thoroughness of structural description are represented in various marked atom strategies for the MA-based descriptors and in toggling the 'dynamic bonds only' option for CGR-based fragments (see section2.1.4). The most detailed structural description was involved herein: correspondingly, the MA3 strategy and local and nonlocal CGR-based fragments. The length of a chosen topology varied from 1 (the fragment standing for atom's count) to 8 for sequences and from 1 to 4 for atom-centered

Figure 36 .

 36 Figure 36. The distribution of the rate constant (log k) for the training set of 8056 Sn1 reactions.

Figure 37

 37 Figure 37. Performed structural encoding of Sn1 reactions by MAbased and CGR-based ISIDA fragment descriptors with the examples of generated fragments of different length (right).

Figure 38 .

 38 Figure 38. GTM class separation for the training set of 8056 Sn1 reactions. The classes are formed in accordance with the bond that breaks while the dissociation.

Figure 39 .

 39 Figure 39. Predicted vs experimental values of the Sn1 rate constant for the external sets of new reaction conditions (test set №1) and new structural transformations (test set №2).

Figure 40 .

 40 Figure 40. Web-implementation of the predictive model for halogen bond strength.

Figure 41 .Figure 42 .

 4142 Figure 41. Web-implementation of the predictive model for tautomeric equilibrium constant.

2 .

 2 SVM model building parameters emerging from the evolutionary libSVM parametrization strategy.

  8.1.1 Data preparation ..................................................................... 8.1.2 Descriptors ........................................................................... 8.1.3 Building and validation of the models ............................................ 8.2 Results and discussions ..................................................................... 8.2.1 GTM visualization of the training set............................................. 8.2.2 Cross validation of the SVR models .............................................. 8.2.3 External validation of the SVR models ........................................... 8.3 Conclusion .................................................................................... 9. Models implementation. ...........................................................................

Conclusion ............................................................................................... References ............................................................................................... Appendix. Part I ................................................................................. Appendix. Part II ................................................................................ Appendix. Part III ...............................................................................

Table 1 . Information about the predictive models developed in this work Studied system Modeled property Encoded information Local descriptors

 1

	Training	Machine-
	set size	learning
		method

a Marked Atoms (MA) based and b Condensed Graph of Reaction (CGR) based local descriptors

Table 2 . The modeling performance of log KBI2 prediction obtained in 5-fold cross validation on the training set of 598 molecules.

 2

Table 3 . Comparison of hydrogen and halogen bond strength. a, b are parameters of linear correlation logKI2=a+b logBHX, s -standard deviation.

 3

  attributed with the experimentally estimated ΔG values for 1:1 complexation. The measurements were carried out in 17 different organic solvents at the temperatures varying from 293K to 303K. From this data, a homogeneous set of 3373 complexes, where the measurements have been carried out in CCl 4 under the standard temperature of 298K, has been extracted and constituted the training set. The set underwent cleaning and filtering excluding all inorganic, metalorganic, deuterium containing compounds and salts. The donors and the acceptors structures underwent a prescribed standardization procedure used on our virtual screening web server (http://infochim.u-strasbg.fr/webserv) and based on ChemAxon's Standardizer 152 (neutralization, isotopes removal, conversion to 'basic' aromatic form etc.). The labeling of the active sites of donor molecules have been performed by SMARTS-based substructure search (by means on an in-house tool using the ChemAxon substructure search API), whereas the acceptors' active sites have been detected and marked by means of the previously developed HB acceptor strength model 113 (http://infochim.ustrasbg.fr/webserv/VSEngine.html). If multiple centers were found, the one with the strongest acceptor propensity was the kept working hypothesis.

Both external test sets have been collected from the same literature. Monofunctional test set (test set 1) consisted of 629 complexes with single HBs. Unlike for the training set, where the strength of the complexes was referred as the Free Energy (ΔG, kJ/mol) for 1:1 complexation at unified conditions, for monofunctional test set the experimental ΔG values

  .5) and on previous experience, showing it to be the best suited in similar contexts113 . The fragmentation schemes for the SVR calculation included various atom coloration by elements symbol, by CVFF force field label or by pharmacophoric types, so that the descriptors are enhanced with an additional chemically relevant information (see section2.1.4). The considered fragments topologies were sequences and atom-centered fragments of the minimal length from 2 to 4 and the maximal length from 3 to 14. Overall, 64 descriptor sets have been tested, where 40 of which, producing the individual SVR models of maximal robustness constituted the consensus SVR model. For MLR, only one atom coloring scheme (by elements) and one fragments topology (sequences) were used. This resulted in 40 descriptors sets used to build 480 MLR individual models, the best of which were kept for the consensus prediction. Cl 4 , C 6 H 6 , C 6 H 5 Cl, CCl 3 CH 3 , C 2 H 4 Cl 2 , C

The prepared consensus SVR and MLR models have been validated on two external test sets.

For the first, monofunctional one, six linear correlations, that relate the ΔG value in a certain solvent and in CCl 4 have been prepared: correspondingly, for C 2 6 H 12

(App., part II, Table II.1-2)

. The Pearson correlation coefficient for the correlations varies from 0.78 to 0.98. Test sets featuring novel combinations of trainingset donor with a training-set acceptor are typically easier to predict than test sets in which a training-set donor is challenged to interact with a never-encountered acceptor or vice versa.

Figure 19. Predictive performance of the consensus SVR and MLR models achieved in 3-fold cross-validation on the training set of 3373 hydrogen bonded complexes.

	DONOUT	87	19	17	4.00	5.09	0.32	0.13
	entire test set 1	629	74	69	3.20	3.81	0.74	0.65
		SVR				MLR
	Class	Number of	Number of outliers	RMSE, kJ/mol	R 2
		compounds	SVR	MLR	SVR	MLR	SVR	M LR
	PAIROUT	262	6	4	2.17	2.51	0.87	0.83
	BOTHOUT	23	1	2	2.18	3.22	0.72	0.41
	ACCOUT	257	48	40	3.83	4.36	0.72	0.64

Table 4 . Predictive performances of consensus SVR and MLR models for different subsets of the external test set №1 possessing single hydrogen bonds.

 4 The donors and acceptors in PAIROUT complexes did occur in the training set but in different combinations, this explains a good models performance for this class. Relatively small RMSE values observed for BOTHOUT subset might be biased by its small size (only 23 H-bond complexes) and composition: thus, most of the compounds in this subset were measured in CCl 4 which decreases the inaccuracies linked to data rescaling from one solvent to another. The largest number of outliers were detected for ACCOUT and DONOUT subsets. Most of wrong predictions in the ACCOUT subset correspond to rare or weak acceptor centers not occurred in the training set, specifically, C-H aromatic acids, compounds with halogen atom acting as acceptor and unsaturated aliphatic or cyclic compounds with double or triple bond acting as acceptor. The worst results obtained for the DONOUT subset could be explained by relatively poor diversity of the donor's class in the training set -over 70% of donors are phenols. For this reason, almost one third of the compounds of this subset are found out of AD.

Figure 20. Predicted (ΔGpred) vs experimental (ΔGexp) free energies for the 1:1 complexes with two cooperative hydrogen bonds. Predicted values were estimated by eq. 37..

  94, RMSE=1.50- 2.22), however still, the MLR method noticeably concede in accuracy. That could be explained, at first, by the lack of information-rich 'colored' descriptor spaces i.e. the ISIDA descriptors enhanced with an additional chemical information (formal charge, force field types, etc.) and, at second, by the important role of non-linearity in HB affinity modeling: thus, out of 40 winning individual SVR models only two of them are based on linear kernel, nevertheless bearing the lowest rank of the statistical score. Successfully cross-validated SVR and MLR consensus models have been challenged for the prediction of two external test sets, the first one of which consisted of complexes with single HBs, measured in either standard CCl 4 , or in other solvent, whereas the second test set was constituted by 12 complexes with cooperative HBs. Apart from the standard model's performance evaluation encompassing all the objects of the test set, a fractional model validation has been performed. The latter considers the test set to be composed of four distinct classes differing by the proportion of 'novelty' of the included compounds. It has been shown that the external sets featuring known partners in novel combinations are indeed easier to predict than sets containing either donors or acceptors that were never seen at the training stage. Logically, the situation should be even more tense for the challenge of predicting a set in which neither acceptors nor donors were met at training stage -however, since that collection was rather small and biased, it was predicted well. The solvent corrections performed with the obtained linear relationships and involved in the assessment of ΔG values of the test set 1 are shown to be a useful tool for the evaluation of ΔG for different solvents, not occurred in the training set. The overall performance referred to the entire test set of single HBs is reasonable: R 2 =0.65-0.74, RMSE=3.20-3.81kJ/mol. At last, on the example of the test set of polyfunctional molecules with multiple intermolecular interactions, it has been shown that the sum of HB affinities for each individual interaction robustly correlates with the observed experimental value (R 2 =0.65-0.87, RMSE=1.63-2.68 kJ/mol), which opens a perspective for the model to be applicable for supramolecular crystal engineering and drug design.

Table 5 . Composition of the DB. The types of tautomerism and the number of transformations in the DB for each tautomeric type.

 5

  ) tautomers not in training set to be predicted under conditions already met among training examples, (3) tautomers in training set to be predicted under novel conditions and eventually, (4) novel tautomers under not yet encountered conditions. Only scenarios (2) and (3) were envisaged here, since not enough external data to support the other two was available. Thus, two test sets, collected from the same literature170, 185 , have been used for external validation. The first test set (test set 1) consists of20 tautomeric transformations (App., part III, Table III.1). which have been occurred in the training set, but under different

reaction conditions. Test set 2 consists of 26 unique transformations (App., part III, Table

III

.2). without structural duplicates in the training set.

Marked atom strategy Number of descriptors Г-score

	MA0	431	0.55
	MA1	204	0.74
	MA2	232	0.73
	MA3	662	0.52

Table 6 . The separation quality of the classification GTM landscapes of the unique structures subset expressed by

 6

Г-score.

Table 7 . The comparison of the performance of SVR and GTM methods for four marked atom strategies. The individual model based on ISIDA atom-centered fragments (length 1-3).

 7

	Marked	Number of		GTM		SVR
	atom	descriptors	Г-score	R 2	RMSE	R 2	RMSE
	strategy
	MA0	445	0.68	0.72	0.82	0.77	0.76
	MA1	218	0.76	0.83	0.64	0.81	0.68
	MA2	246	0.75	0.84	0.63	0.82	0.67
	MA3	676	0.69	0.78	0.73	0.80	0.71

Table 8 . Example of tautomeric transformations with considerable log KT difference for different solvents.

 8

  (4+2) cycloaddition, about 8% (3+2) dipolar cyclisation, and 7% (2+2) cycloadditions. The measurements were carried out in C 6 H 5 CH 3 , CH 3 COOC 2 H 5 , C 6 H 5 OCH 3 , C 6 H 5 NO 2 , CH 3 CN, C 2 H 5 OH, C 6 H 5 Cl, THF, CH 2 Cl 2 , DMSO, CHCl 3 , C 2 H 4 OC 2 H 4 O, C 5 H 11 OH, C 6 H 5 Br, C 6 H 6 , CH 3 COCH 3 , C 6 H 12 and CH 2 ClCH 2Cl solvents at temperatures varying from 273 to 423 K. The most frequently occurred dienes for the Diels-Alder series are: condensed aromatic cycles, cyclopentadienones, cyclopentadienes and benzofuranes. Among dienophiles, the prevailing structures are: maleic acid derivatives, cyanoethylenes, Phsubstituted ethylenes and benzoquinones. The data set (reagents and products) underwent an accustomed standardization protocol using ChemAxon's Standardizer Utility: basic aromatization, isotopes removal, NO 2 -, NO-, N 3 -, RRSO 2 -, CN-group transformation152 While cyclization due to the reaction mechanism different stereoisomers can be formed. An analysis of the rate constants of reactions forming two different stereoisomers have shown that logk difference was from 0.01 to 1.0 log values. Based on our previous work209 we considered that the difference is too small to be taken into account and comparable with the interlaboratory errors. Thus, for the reactions able to form different stereoisomers during cyclisation, the logk constants were calculated as a mean of given experimental values. No other duplicates were found in the dataset. The 1849 reactions remaining after cleaning have been characterized by 1849 values of logk, 1356 values of Ea (kJ/mol) and 1237 values of logA, where among the latter there were no reactions without Ea. The prepared data has been divided into the training and test sets. External test set. The test set consisted of 200 reactions randomly picked from the initial data set. For the reactions possessing logA and Ea values and comprised into the test set, the corresponding values were not taken into consideration. In such a way, the test set is considered to be attributed with 200 experimentally measured logk values solely. Out of 200 reactions in the test set, there were 57 structurally new transformations, that did not occur in the training set. Individual statistics for them is discussed (see

Table 9 . Examples of reactions projected into the zones corresponding to Figure 29.

 9

Table 10 . Predictive performance of the SVR models, built separately for each of the property, and a single GTM model, universal for all three properties, in 5-fold cross validation for the training set of 1649 reactions of cycloaddition. The descriptor space is based on atom-centered fragments of the length 2-3.

 10

Table 11 . Validation of four different methods of logk assessment on the test set of 200 cycloaddition reactions. Notice that only 164 reactions are retained by the Fragment Control applicability domain.

 11

  Since the external test set has been randomly chosen from the initial data set, it includes reactions already occurred in the training set but proceeding under different reaction conditions. That may lead to an overestimation of model's predictive performance, due to occurrence of same reactions in both training and test sets. Thus, individual statistics for the structures, never encountered in the training set provides an unbiased assessment of the model's predictive ability with regard to new structures. Among 200 external set's reactions, only 57 transformations didn't occur in the training set. Table11shows that only direct logk assessment with SVM model leads to reasonable correlation between predicted and experimental values: R 2 = 0.72 which is close to the determination coefficient on the validation stage for the entire test set. On the other hand, application of the Arrhenius-based SVR assessment as well as GTM-based models result in relatively poor predictions R 2 < 0.5.

	Subset of unique structural transformations.
	Log k predictive	SVR		GTM
	method	R 2	RMSE	R 2	RMSE
	Direct	0.72	0.80	0.38	1.21
	Arrhenius-based	0.33	1.26	-0.48	1.87
	Temperature-	-	-	0.40	1.15
	corrected GTM
	GTM-based direct logk		GTM-based temperature-
	prediction		scaled approach
		R 2 =0.90			R 2 =0.94
		RMSE=0.48		RMSE=0.39

Table 12 . The comparison of the direct , Arrhenius-based and temperature-scaled logk assessment for the subset of 57 unique structural transformations. Figure 33. Performances of SVR and GTM models for the subset of unique structural transformations (57 reactions) obtained with the direct log k prediction method. The statistics is given in Table 12.

 12

  The measurements have been carried out in: CH 2 OHCHOHCH 2 OH, CH 2 ClCH 2 Cl, C 2 H 4 OC 2 H 4 O, CF 3 CH 2 OH, CH 2 ClCH 2 OH, (CH₃)₃COH, CH 3 COOH, CH 3 CN, C 2 H 5 COHCH 3 , C 4 H 9 OH, CH 2 OHCH 2 OH, C 2 H 5 OH, C 2 H 5 OC 2 H 5 , NH 2 COH, D 2 O, C 6 H 11 OH, DMSO, CH 3 OH, C 6 H 5 NO 2 , DMF, C 8 H 17 OH, THF, C 3 H 7 OH, CH 3 CHOHCH 3 , CH 3 COCH 3 , H 2 0, C 4 H 8 SO 2 and C 6 H 5 CH 3 . An accustomed standardization protocol of ChemAxon's Standardizer Utility has been applied: basic aromatization, NO 2 -, NO-, N 3 -, RRSO 2 -, CN-group transformation 152

Table 13 . Predictive performance of the best individual SVR models in 5-fold cross-validation for the training set of 8056 SN1 reactions.

 13

Table I . 3 .

 I3 Predictive performances of the MLR consensus models in 5-fold cross-validationinvolving the different marked atom strategies using bounding box and fragment control AD approaches.

	Mark atom	MLR CM with AD
	strategy	5CV RMSE	Q 2
	MA0	0.46	0.902
	MA1	0.39	0.930
	MA2	0.36	0.944
	MA3	0.36	0.942

Table I . 4 .

 I4 The list of outliers. The third column contains name, solvent and log KI2 value correspondingly.

	Table I.5. Halogen bond versus Hydrogen bond
	Here, logKBI2 is logarithm of stability constant of the
	1:1 complexation of organic Lewis bases with I2 in
	heptane, cyclohexane, hexane or methylcyclohexane
	at 298K,
	logKBHX is logarithm of stability constant of the 1:1
	hydrogen bonding of the same organic bases with 4-
	F-phenol in CCl4 at 298K.
	Structure	Molecule Title
	1	Tetramethylguanidine, Hept, 4,37
	2	N,N-Diisobutylmethylamine, Hept, 1,45
	3	2-Dimethylaminopyridine, cHex, 0,95

Table II . 2 .

 II2 Linear Gibbs energy relationships for the Gibbs energy of H-bonding: the parameters of the linear equation G (in CCl4) = a G (in solvent) + b, where R is Pearson correlation coefficient, the temperature range is 293 -298 K

	no.	solvent	a ± a	b ± b	R 2	s	n
	1	C2Cl4	0.975 ±0.027	0.42 ±0.29	0.982	0.70	25
	2	C6H6	1.024 ±0.034	1.44 ±0.30	0.916	1.62	87
	3	C6H5Cl	1.089 ±0.062	0.02 ±0.66	0.951	1.13	18
	4	CCl3CH3	0.962 ±0.066	2.14 ±0.89	0.922	1.36	20
	5	1,2-C2H4Cl2	1.109 ±0.063	3.07 ±0.51	0.783	1.95	86
	6	C6H12	0.891 ±0.020	0.64 ±0.19	0.929	1.71	145

Table II . 3 .

 II3 Experimental and predicted by SVM model complexation free energies ∆G of the external test set with two cooperative H-bonds (dimers set).

Table II . 4 .

 II4 SVM model building parameters emerging from the evolutionary libSVM parameterization strategy.

	Descriptor space	LibSVM setup	3CV RMSE	3CV Q 2

Table II . 5 .

 II5 Predictive performances of the MLR consensus models in 3-fold cross-validation involving the different marked atom strategies without accounting for the models applicability domain.

	Descriptor	MLR CM
	strategy	3CV RMSE	3CV R 2 det
	MA0	2.73	0.796
	MA2	2.48	0.832
	MA3	2.30	0.855
	Here R 2	det is the squared determination coefficient of 3CV predictions

Table II . 6 .

 II6 Predictive performances of the MLR consensus models in 3-fold cross-validation using bounding box and fragment control AD approaches. The MA3 type of descriptors were used. the number of species within D. R 2 det is the squared determination coefficient of 3CV predictions

		MLR CM with AD
	3CV RMSE	3CV R 2	det	npred a
	2.11	0.880		3084

a

Table II . 7 .

 II7 Predictive performances of the MLR Consensus Model for different classes of the external test set accounting for the model's applicability domain.

	Class

name Number of compounds Number of predicted compounds RMSE R 2 det

  experimental G values were recalculated to CCl4 using linear correlations (see TableS1)

	10 16				20 50		Water Acetone	0.13 -0.12/-0.01 0.01 -0.27/-0.32
								(97%)
	70				-25.86		-14.07		CCl4	ACCOUT
	11 17				20 20		Ethanol Water	-1.53 -1.34/-0.21 -0.39/-0.34
								(61%)	0.92
	71				-4.18		1.25		CCl4	DONOUT
	12 18				30 49.8		Methanol DMSO	-0.6 -	-0.9/-0.77 -0.01/-0.31
										0.31
	72 13 29			-6.26	30 25.1	2.49 1-Propanol C6H12 -1 Tetracloro -	DONOUT -1.01/-1.27 -0.17/-0.07
								methane	0.13
	14 20				10 37		Toluene Tetracloro	--	-0.09/-0.22 -0.01/-0.19
	73				-13.81		-7.75 methane	CCl4 0.32 0.69	PAIROUT
	PAIROUT 15	262	262		2.51 40	0.830 Toluene	0.96 -0.12/-0.01
	BOTHOUT	23	11		3.63		-
	74	ACCOUT	257	105	-13.30	3.34 -3.57	0.735	CCl4	DONOUT
		DONOUT	99	44		4.35		0.164
	Solvent CCl4	287	179		3.63		0.637
	Solvent Other	354	243		2.45		0.843
	entire test set	629	422		3.01		0.759

a b SVM method has been used. c see annotations in section 2.1.2.

Table III . 2 .

 III2 External test set 2.

	1 14		20 20 60	Water Ethanol BromoB	3.48 3.46/3.18 out 0.6 -in 0.27 0.18/0.06 in
				enzene		0.10/0.01
	2		26	Acetone	0.31 0.14/0.86 out
	15		20 35 Acetonitr Acetone	0.19 -0.17	0.65/--0.01/-	out in
				ile		0.08 0.56
	3		35	Water	0	-1.47/-	out
	16		20 35	Dioxan Water 1,4-	0.2 0.96 -1.39/-1.42 0.39/0.33 in -1.67 in
	4		60 Nitroben	1.38	1.19/1.2	in
				zene
	17		25 70	Ethanol DMSO	0.03 -1.28 0.51/0.47 in 0.13/-in
						0.03
	5		20	Water	3.7	-2.39/-	out
	18		10 100	Toluene -0.22 Water -6.96	0.44 -0.21/--2.19/-	in out
						0.01 2.55
	6		20	Water	0.69	-1.34/-	in
			100	Water	-6.95	0.2 -2.19/-	out
	19		25	Water	0.31 0.34/0.21 in 2.55
	№ 7	External test set №2	T, °C 35 100	Solvent Water Water	Log KT exp 1.99 -0.69	Log KT SVM/GT 1.91/0.25 0.27 pred --0.88/-	AD out out
						M

Table III . 5 .

 III5 SVR individual models constituting the web-deployed consensus model. Model building parameters emerged from the evolutionary libSVM parametrization strategy.

Thus, all the considered types of local descriptors could be used for structural representation of chemical processes, however, some of them have distinct shortcomings. Thus, the substituent constants, coming from the experiment and measured only for some of the groups, could not be served as a universal description and suitable only for homogeneous data sets of few varying substituents. The quantum-chemical descriptors, correspondingly, require high computational costs, which is a strong limitation in case of modeling of big data sets. The electrotopological indices are not expensive computationally, however, they do not reflect the structural aspects explicitly, and are composed for each atom as the sum of the corresponding electrotolopogical aspects, thus becoming not interpretable in terms of chemical structure. Moreover, neither of these descriptor types supports an explicit emphasis on the active centers. Therefore, ISIDA descriptors offer the best solutions to the abovementioned constraints due to fast computation, direct structural representation and various possibilities of designation of the active sites.

Acknowledgements

I would like to express my great appreciation to my supervisors, Prof. Alexandre Varnek and Dr. Timur Madzhidov for their thorough guidance, continuous solicitude and care. Besides of my supervisors, I am particularly grateful for the assistance given by Dragos Horvath, who directed me in the practicalities of the field, carefully taught and instruct me all over my PhD studentship. I also wish to acknowledge the help and advises provided by Gilles Marcou. I

The definition of the AD is a crucial aspect since the prediction of an object which is being outside the AD is unreliable and could lead to wrong conclusions and undesired consequences.

A lot of different schemes proposed for the AD determination, that could rely purely on the descriptors constituents or could be derived from the machine learning method. The designation of the most appropriate AD is still a matter of discussion.

The AD of all the projects of this study is based on Bounding Box, for each descriptor vector reckoning the minimum and maximum values encountered in the training set. An object is considered to be out of AD if at least one of its descriptor values violates the defined min-max range. 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 13. Representation of the concept of applicability domain for the chemical space based on two descriptors. The prediction of the test compound inside the domain is reliable whereas of the compound outside the domain the prediction is not trustworthy.

Chapter 6

QSPR modeling and visualization of tautomeric equilibria.

The complex phenomenon of tautomerism is still a challenge for chemoinformatics and computational chemistry in terms of quantitative estimation, mode of transformation and representation. Tautomerism is ubiquitous and plays a key role in practically important processes including biochemical ones, such as the relation of tautomeric transformations to spontaneous mutations as a consequence of mispairing by rare tautomeric forms of purines and pyrimidines [171][172][173][174] and its relation to enzyme-substrate interactions [175][176] . As for organic chemistry, the prevailing conformation of a certain tautomeric form may affect the product in a chemical reaction. The field of drug design likewise needs the determination of the predominant ligand structure for virtual screening and modeling [177][178] . The number of works devoted to the elucidation of the qualitative and quantitative aspects of tautomerism is still being extended [178][179][180][181] . However, most of them are related to quantum chemistry. In spite of the importance of this phenomenon, only two software tools dedicated to the assessment of the tautomeric population are currently available: the Marvin Tautomerization Plugin 182 and TauThor/MOKA 183 . Both of those tools estimate the equilibrium constants in water at room temperature using predicted pKa values for all individual tautomeric forms. In many important cases, their predictive performance appears, however, to be too low, because of Chapter 9

Models implementation.

Halogen bond basicity of organic molecules.

The model predicting the strength of halogen bonding between an organic molecule and diiodine (log K I2 ) is available on the web-server http://infochim.ustrasbg.fr/webserv/VSEngine.html. The SVR consensus model consist of 9 best individual models with fitness score (R 2 ) from 0.903 to 0.920. The detailed information about the employed descriptors spaces is given in Appendix (part I, Table I.2).

The model allows:

• An automated detection of the main halogen bond acceptor