Arkadii Lin

Natalia Dyubankova

Timur I Madzhidov

Ramil I Nugmanov

Jonas Verhoeven

Timur R Gimadiev

Valentina A Afonina

Zarina Ibragimova

Asima Rakhimbekova

Pavel Sidorov

Andrei Gedich

Rail Suleymanov

Ravil Mukhametgaleev

Joerg Wegner

Hugo Ceulemans

Alexandre Varnek

Atom-to-atom mapping: a benchmarking study of popular mapping algorithms and consensus strategies

Keywords: chemical reactions, data cleaning, atom-to-atom mapping

 which implement different AAM algorithms. An open-source RDTool program was optimized, and its modified version ("new RDTool") was considered together with several consensus mapping strategies. The Condensed Graph of Reaction approach was used to calculate chemical distances and develop the "AAM fixer" algorithm for an automatized correction of erroneous mapping. The benchmarking calculations were performed on a Golden dataset containing 1851 manually mapped and curated reactions.

The best performing RXNMapper program together with the AMM Fixer was applied to map the USPTO database. The Golden dataset, mapped USPTO and optimized RDTool are available in the GitHub repository https://github.com/Laboratoire-deChemoinformatique.

Introduction

Atom-to-atom mapping (AAM) [6,7] is a procedure that establishes a correspondence between the atoms of reactants and products. AAM allows to identify a reaction centre (RC) which, in turn, helps to prepare reaction templates used in an automatized forward/retrosynthesis planning [8][9][10][11][12] , as well as to perform reaction classification [START_REF] Hendrickson | Encycl. Comput. Chem[END_REF] and reaction searching [START_REF] Lin | [END_REF]15] .

Several publicly and commercially available AAM tools are currently available. Their comparison would allow choosing an appropriate tool suitable for a given project. Some benchmarking studies have been already reported, [16,17] but they became outdated once new mappers were released.

In this paper, we compare the most popular AAM tools: Automapper (ChemAxon) [START_REF]can be found under[END_REF] , Indigo [START_REF] Savelev | [END_REF] , RDTool [START_REF] Rahman | [END_REF] , NameRXN (NextMove) [START_REF] Software | NameRXN," can be found under www[END_REF] , and RXNMapper [START_REF] Schwaller | [END_REF] which implement different AAM algorithms. Thus, the opensource Indigo toolkit [START_REF] Savelev | [END_REF] applies the VF2 algorithm [18] which is an advanced graph matching algorithm based on Maximum Common Substructure (MCS) isomorphism principle [19] . This algorithm detects common structural motifs in reactants and products followed by identification of related pairs of atoms.

ChemAxon's Automapper uses a combination of MCS and a minimal Chemical Distance (CD) approaches. According to Chen et al. [20] , a chemical distance is calculated as a sum of two terms: the number of chemical bonds created, broken, or transformed, and the number of atoms, which charge or valence state has been modified in chemical reaction. The AAM procedure with ChemAxon's Automapper consists of two steps: 1) detection of an appropriate MCS, and 2) CD-based optimization of mapping the atoms that are not a part of the MCS.

The NameRxn software [START_REF] Software | NameRXN," can be found under www[END_REF] from NextMove performs AAM using large ensemble of rules describing known reaction mechanisms/transformations. The algorithm classifies a given reaction according to the implemented heuristics, followed by identification of the reaction centre described by an ensemble of atoms at chemical bonds formed, broken or transformed in course of chemical RXNMapper tool [START_REF] Schwaller | [END_REF] implements a deep learning algorithm that is much faster than the abovementioned ones because it does not need to perform time-consuming combinatorial-seeking operations. RXNMapper uses a model based on the transformer architecture [21][22][23] which relates SMILES strings [24] of reactants (input) and products (output). An attention mechanism is used to assign a correspondence between the atoms of reactants and products. RXNMapper was trained in unsupervised settings on a subset of the USPTO database [25] for which the AAM accuracy >99% was reported.

Reaction Decoder Tool (or RDTool) uses a consensus approach combining four distinct MCS-based algorithms: Mixture-MCS model, Min-sub model, Max-sub model, and Assimilation model [26] . All of them are based on pairwise MCS identification between reactant and product molecules. In the Min-sub model (further mentioned as MIN), the smallest MCS between the reactants and products are first mapped. On the contrary, in the Maxsub model (further mentioned as MAX), the largest MCS between the reactants and products are first mapped. Then, the mapped part is removed from the reaction followed by the repetition of pairwise MCS identification until no more atoms can be mapped. The user-defined structural matches instead of MCS can be applied in preference. The Mixture-MCS model (further mentioned as MIXTURE) shares some features from MIN and MAX approaches. Assimilation model uses a restricted MCS algorithm that prefers the ring-to-ring mapping over the ring-to-chain one. Among those, the best AAM solution is selected according to a score accounting for (i) fragments count, (ii) dynamic bonds count, and (iii) energy of bonds cleavage/formation. According to Rahman et al. [START_REF] Rahman | [END_REF] , the best solution corresponds to the minimal structural changes. In line with this principle, the algorithm minimizes the bond changes in the reaction and, therefore, reduces the chemical distance and energy costs for forming or breaking the bonds. To compute the energy of bonds cleavage/formation, the precomputed values are used [START_REF] Stevenson | The Strengths of Chemical Bonds[END_REF][START_REF] Deb | [END_REF][START_REF] Benson | [END_REF] . Notice that the number of dynamic bonds is counted using the Dugundji-Ugi matrix [START_REF] Dugundji | Comput. Chem[END_REF] , which requires a fully balanced reaction equation. For this reason, RDTool performs some data pre-processing including reaction balancing. This, unfortunately, cancels previously performed reaction standardization.

It should be noted that any mapper assigns to the atoms the numbers (hereafter, called "AAM numbers") further used in the reaction SMILES or RXN file to identify atom-to-atom mapping. Among the benchmarked mappers, only ChemAxon mapper uses the initial ordering of atoms as AAM numbers.

All benchmarking calculations were performed on a reference Golden dataset containing 1851 manually mapped and curated reactions.

Method

Optimized RDTool and consensus solutions

Since RDTool is an open-source software, we decided to refactor the code in order to use independently three out of four constituent algorithms: Min-sub (MIN), Max-sub (MAX) and Mixture-MCS (MIXTURE) models. The fourth algorithm was omitted due to the code fuzziness. The code was also revised in order to accelerate the functioning of the program. This new version of the program was named new RDTool.

In complement to individual AAM tools, six consensus scenarios were considered:

1

)

Assessment of AAM accuracy with Condensed Graph of Reaction

Chemical Distance corresponding to a given AAM solution can be easily assessed for any reaction transformed into a Condensed Graph of a Reaction (CGR) [START_REF] Hoonakker | [END_REF]32] . CGR is a single graph, which encodes an ensemble of reactants and products as shown in Figure 1. CGR results from the superposition of the atoms of products and reactants having the same numbers. The nodes in CGR correspond to "static" or "dynamic" atoms whereas the edges correspond to "static" or "dynamic" bonds. The state of dynamic objects -charge and radical label for atoms and bond order for bonds -is changed in course of reaction. Once the AAM operation is performed, a related CGR can be prepared even if a reaction equation is not balanced (e.g., a minor reactant or a leaving group is missed). In line with the Chen et al. [20] definition, a Chemical Distance for a given reaction is estimated by the number of related dynamic objects. Clearly, different AAM solutions lead to different CGRs. In this study, CGRs have been prepared for all reactions from the reference set according to AAM performed either manually by experts in organic chemistry or using considered software tools. Ideally, a mapping provided by a particular tool corresponds to the manual one, i.e., the related CGRs should be identical. Notice that in some cases, erroneous mapping leads to the same chemical distance as that issued from the manual AAM (see example in Section 2.3). Therefore, the identity of two condensed graphs is a more solid criterion of the AAM correctness than a simple comparison of related chemical distances.

Heuristics for correction of erroneous AAM

Systematic errors of atom-to-atom mapping can be corrected using heuristics based on the CGR approach. The idea is to replace a structural motif representing an erroneous reaction centre and its environment in a CGR with another one corresponding to a correct mapping using empirical remapping rules. The concept of AAM fixing was proposed in our early publication [33] ; here we demonstrate its fully automatized implementation.

The overall process can be split into two parts: (i) fixing rules preparation (see Figure 2), and (ii) fixing rules application (see Figure 3). A fixing rule is generated for a pair of correct/erroneous AAM generated for a given reaction and represented in either RXN or SMILES formats. A fixing rule is generated in 4 steps:

1. Two related CGRs are prepared for correctly and erroneously mapped reactions. For each CGR, a reaction centre with its closest environment (RCclose) is extracted (see Figure 2a). 2. Preparation of "extended reaction centres", RCext, by augmentation of correct RCclose fragment by the atoms belonging to the erroneous RCclose and vice versa. In such a way, both correct and erroneous RCext fragments contain the same set of atoms (Figure 2b).

3. In order to determine the required permutations P of AAM numbers, the extended reaction centres are deconvoluted to reaction equations.

A permutation rule can be obtained by overlapping the atomic AAM numbers' lists extracted from the corresponding deconvoluted reaction equations (Figure 2c). 4. A fixing rule is compiled out of three instances:

an erroneous and correct RCext and a permutation rule (Figure 2d).

Notice that the order of atoms in the two atoms' lists at step 3 must be the same: the algorithm crashes or generates wrong permutation solutions when the atoms are shuffled by the mapper. Among all benchmarked here tools, only ChemAxon mapper keeps the order of atoms as in the initial RXN of SMILES. Therefore, AAM Fixer was applied only to the reactions mapped by ChemAxon.

Once a fixing rule is prepared, it can be applied to fix AAM numbers for a given reaction for which a CGR is generated (Figure 3a). Then erroneous extended reaction centre is embedded in the CGR (Figure 3b). Hits found for a particular RCext signal that the AAM numbers are wrong and can be fixed with the corresponding permutation rule P (Figure 3c). AAM numbers in the product part of the reaction are permuted according to P. Finally, a new CGR is generated for the corrected reaction followed by the embedding of the correct RCext in order to check whether the AAM fixing was done successfully (Figure 3d). The original AAM numbers are restored if the latter embedding fails. Once extracted, the fixing rules can be applied to correct AAM generated by any mappers.

Fixing rule extracted from a given reaction is validated by its application to the same reaction. In some cases, the correct mapping cannot be achieved. Usually, it concerns the cases when atoms distant from the reaction center need to be remapped, see Figure S1 in SI. Finally, all validated fixing rules are integrated into AAM Fixer. (d) a new CGR is generated, and the correct RCext is embedded to check if the AAM correction is done successfully. Here, the dynamic bond "1>0" corresponds to a broken single bond, and "0>1" corresponds to a formed single bond.

Data

The data set used in this study as a reference resulted from the merging of curated dataset used in the benchmarking study by Jaworski et al. [34] and a subset of reactions extracted from USPTO database. The first set contains reactions with full stoichiometry and taken from the Organic Syntheses collection [START_REF]Organic Syntheses[END_REF] , randomly selected and typically mostly stoichiometrically unbalanced reactions from Reaxys [START_REF] Goodman | [END_REF] , mechanistically complex reactions taken from various literature sources (e.g., Kurti's "Strategic Application of Named Reactions in Organic Synthesis" [START_REF] Kurti | Strategic Applications of Named Reactions in Organic Synthesis[END_REF] or Grossman's "The Art of Writing Reasonable Organic Reaction Mechanisms" [START_REF] Grossman | The Art of Writing Reasonable Organic Reaction Mechanisms[END_REF] ), relatively simple reactions with full stoichiometry taken from total syntheses published in Org. Lett., J. Am. Chem. Soc., and J. Org. Chem., typical reactions without full stoichiometry extracted from patents, and mechanistically complex reactions which include rearrangements and multicomponent reactions taken from recent literature (in most cases, after 2010 and from Org. Lett., J. Am. Chem. Soc., and J. Org. Chem.). A detailed description of reactions types included into the dataset can be found in ref. [34] The second set extracted from USPTO gathers examples of the most popular reaction types used in medicinal chemistry (e.g. SN2, Buchwald-Hartwig, Suzuki-Miyaura coupling, etc.) [START_REF] Brown | [END_REF] .

A manual check of 1405 organic reactions from the Jaworski's set [34] revealed several serious problems with structures preparation or atom-to-atom mapping. Thus, 10 reactions were identified as duplicated (see an example in Figure S2 in Supporting Information), 30 reactions were discarded since they contained molecules with invalid valences (Figure S3) or they could not be interpreted (Figure S4). The discarded reactions were analyzed manually and curated, if possible. Namely, 17 reactions with chemical structure representation errors were corrected. We've also discovered some 350 reactions containing radicals that were manually curated (see an example of such reaction in Figure S5).

Manual inspection of atom-to-atom also revealed several problems. Eight reactions containing obvious errors (e.g., one carbon atom in a reactant corresponded to both oxygen and carbon atoms in a product simultaneously) were manually curated (see an example in Figure S6). On the other hand, 22 reactions with nonobvious mapping which may depend on experimental conditions (see Figure S7), and 3 reactions with unclear reaction mechanisms or lost reactants (see Figure S8) were removed. Curated Jaworski's dataset containing 1382 reactions was extended with 469 manually mapped USPTO reactions. The final reference ("Golden") dataset comprises 1851 reactions. Its composition has been analysed with the help of the NameRxn (NextMove) tool. The tool failed to classify 628 reactions. Occurrences of particular reaction types for remaining 1164 reactions are given in Supporting Information (see Figure S9 in SI).

In addition, the entire USPTO database [25] containing some 1.3M unique and standardized reactions [40] (functional groups standardization, dearomatization, valence checking, selection of major tautomers, removal of explicit hydrogens, cleaning radicals and isotopes, clear stereo, split of ions, aromatization, duplicates removal) was used in order to demonstrate the functioning of AAM Fixer.

Results and Discussion

Benchmarking results

15 mapping strategies described in the Method section were applied to the Golden dataset. As it follows from Figure 4, ChemAxon is the fastest tool: it takes approximately 0.003 second to map one reaction. NameRXN (NextMove) and RXNMapper follow ChemAxon with 0.02 and 0.05 seconds per reaction, respectively. RDTool is the slowest one: it takes around 2.46 seconds to map one reaction. This is not surprizing since RDTool executes four mapping algorithms followed by calculation of seven metrics used to perform consensus prediction. Revision of its code significantly accelerated the mapping procedure: the new RDTool including three AAM algorithms (MIN, MAX, and MIXTURE) is about 40% faster than the original tool. Figure 5 and Table 1 report an information about AAM performance of different tools and algorithms assesed for the Golden set. As explained in Section 2.2, AAM performed for a given reaction was considered correct if related CGR coinсided with that obtained by the manual mapping. Compared to other tools, NameRXN [41] demonstrates the lowest rate of correct atom-to-atom mapping (33.6%). This can be explained by the fact that the AAM algorithm uses a limited set of reaction classification rules. If there is no rule characterizing a given reaction, no mapping is performed. In addition, existing rules do not always guarantee a correct AAM (see Figure 6). The RXNMapper tool demonstrates a good performance: 1550 reactions (83.74%) were correctly mapped which is the best result among all considered reaction mappers. Unlike other studied mappers, RXNMapper consistently generates correct AAM for esterification and esters hydrolysis reactions. However, this tool has a serious drawback: it resets any previously performed standardization operation and reorders the reactants and the products (Figure S10).

The original RDTool succeeded to map 76.23% of reactions in the reference data set, which is very similar to the new RDTool (76.18%). The code revision performed in this work allowed us to significantly reduce the time of computations (see Figure 4) and to apply MIN, MAX, and MIXTURE algorithms separately (see Method section). Comparing performances of the new RDTool and its constituting individual mapping algorithms, we found that the MAX algorithm already provides a 73.96% accuracy, which corresponds to 97.1% of the total number of reactions correctly mapped by the new RDTool. Thus, the MAX algorithm is almost as accurate as the entire RDTool but is 5 times faster.

The accuracy of ChemAxon's AAM (70.45%) is similar to earlier reported results. [START_REF] Schwaller | [END_REF]34] Overall, 547 reactions were erroneously mapped.

The success rate of the Indigo mapper is relatively low, especially for halogenation reactions. In some cases, the tool maps a halogen atom in a product as it comes from an unspecified reactant, even if a Hal2 molecule is present in the reactants part (see an example in Figure S11 in SI). In addition, in the case of symmetric molecules, Indigo makes one-to-many correspondence of one reactant atom to several product atoms, i.e., it assigns the same numbers to the atoms in different symmetrically identical fragments (see an example in Figure S12). This leads to significant problems with further usage of mapping results, e.g., with CGR generation. Therefore, reactions with oneto-many correspondence were classified as "unsuitable" (see Figure 5 in red).

Among all examined consensus solutions, only Consensus 3 combining RXNMapper, MAX, MIXTURE, MIN and ChemAxon methods performs slightly better than RXNMapper individually (84.06% vs 83.74%). This negligible accuracy improvement does not compensate a 30 folds slowdown of the consensus compared to the individual solution.

Besides that, Table 1 demonstrates that the lowest Chemical Distance (CD) does not always guarantee correct AAM assignment. Thus, for a large portion of reactions, erroneous and correct mapping correspond to equal chemical distances whatever the mapping algorithm is. Moreover, in some Diels-Alder reactions, the lowest CD corresponds to erroneous AAM (see Figure 7).

Another question concerns the difference in the AAM algorithms implemented in the benchmarked tools. Table 1 shows that the number of reactions correctly mapped with only one selected algorithm is negligible except of that for RXNMapper which is the only tool performing correct mapping for esterification and esters hydrolysis reactions on the Golden dataset. The latter, however, wasn't reproduced for the entire set of USPTO reactions. It is also surprising that according to our tests, the tool perfectly maps reactions of some particular types in the Golden dataset, but it fails on those from USPTO [25] . In particular, this concerns reduction of carboxyl groups to alcohols (89% reactions of this type in USPTO were erroneously mapped), methyl ester formation in methanol (31%), and esterification reactions (10-30%). c Chemical Distance for the erroneously mapped reaction is equal to that with correct AAM.

d See description of the consensus strategies in Section 2.1.

e not relevant for consensus strategies

Correction of erroneous mapping with the AAM Fixer

In order to investigate the efficiency of the AAM Fixer tool, it has been applied to the mapped Golden dataset aiming (i) to extract the fixing rules and (ii) to correct erroneously mapped reactions.

For 547 (out of 1851) reactions incorrectly mapped by ChemAxon, 535 unique fixing rules were extracted. Out of those, only 442 rules have been validated (see section 2.3) which allowed us to correct 452 reactions with erroneous mappings. Interestingly, that the rules derived for ChemAxon mapping allowed us to correct 62 reactions from the Golden set erroneously mapped by RXNMapper. This proves the transferability of the fixing rules. Unfortunately, it is impossible for the moment to extract fixing rules for reactions mapped by RXNMapper since all considered here AAM algorithms except ChemAxon mapper change the atoms' order in the RXN file (see Section 2.3).

Finally, AAM Fixer powered by the rules, extracted for the Golden dataset, was applied to the USPTO dataset mapped by RXNmapper. This led to AAM corrections for 2021 reactions. Such a small portion of fixed reactions is caused by the fact that rules are quite selective due to the RC environment of radius 1 (RCclose).

The advantage of the AAM Fixer is its ability to fix AAM in cases when the chemical distance for a wrong AAM is the same or even lower than for a correct one, like Diels-Alder or esterification reactions (Figure S13 in Supporting Materials). Notice that the fixing rules can be prepared manually and then implemented in the AAM fixer. Here, the reaction itself is on the left side, and the corresponding CGR is on the right side. The symbol "n>m" corresponds to the transformation of bond order n to m [33] . Notice that the chemical distance corresponding to erroneous mapping (CD=5) is lower than that for correctly mapped reaction (CD=6).

Conclusions

Benchmarking study comparing the performance of Atom-to-Atom Mapping of 5 popular tools (RXNMapper, ChemAxon, NameRXN, Indigo, and RDTools) and several algorithms and consensus strategies applied to the "Golden" dataset revealed that no one provides with an ideal AAM assignment. RXNMapper is ranked the best with the accuracy 83.74% and a reasonable speed of calculations (0.05 seconds). However, the deep learning model integrated into this tool looks overfitted because it fails to map certain types of USPTO reactions although very similar reactions from the Golden set were perfectly mapped.

It has been demonstrated that the lowest Chemical Distance doesn't always guarantee correct AAM assignment. Thus, for some Diels-Alder reactions, erroneous mapping corresponds to smaller CDs compared to correct AAM. For a large portion of reactions, erroneous and correct mapping correspond to equal chemical distances for any studied tool. For such cases, we have developed a rule-based AAM Fixer tool that corrects frequent types of AAM errors. It extracts rules from a pair of correct and wrong AAM and can fix mapping for new reactions.

Finally, we've proven that the Condensed Graph of Reaction approach is particularly useful for chemical distances assessment, detection and correction of wrong mapping.

Data availability.

The manually mapped "Golden" dataset containing 1851 reactions, the USPTO database mapped with RXNMapper, standardization workflow and the optimized RDTool (new RDTool) are available on GitHub: https://github.com/Laboratoire-de-Chemoinformatique/.The AAM Fixer tool with embedded rules is available on GitHub: https://github.com/cimmkzn/AAMFixer .

Funding

This project was supported by the VLAIO HBC.2018.2287 MaDeSMArt (Machine Design of Small Molecules by AI) grant. The AAMFixer approach development was supported by the subsidy allocated to Kazan Federal University for the state assignment in the sphere of scientific activities.

Conflict of interests.

There is no conflict of interest

Figure 1 .

 1 Figure 1. The reaction between malonic ester and 4-chloro-beta-nitrostyrene (left) and related Condensed Graph of a Reaction (CGR) (right). Dynamic bonds in CGR are depicted by dash lines. The dynamic bond "2>1" between two carbon atoms 8 and 9 corresponds to the transformation of a double bond to a single bond, whereas "0>1" corr esponds to the formed single bond between atoms 8 and 17. The reaction centre (RC) includes atoms 8, 9, and 17 at the dynamic bonds.

Figure 2 .

 2 Figure 2. AAM fixing rule preparation workflow illustrated on the example of the addition-elimination reaction for which erroneous mapping was obtained with ChemAxon. (a) The RCclose fragments corresponding to the correct and erroneous mapping include six atoms only four of which (atoms 2, 3, 10 and 13) are in common. (b). Eight atoms, that belong to correct or wrong reaction centres or shared by both (atoms 1, 2, 3, 4, 9, 10, 11 and 13) are considered in RCext. Thus, extended reaction centres are prepared by the addition of atoms 9 and 11 to the correct RCclose fragment and of atoms 1 and 4 to the erroneous RCclose. (c) Deconvolution of RCext to reaction equations and superposition of the related atoms in the products reveals a permutation rule. ( d) The numbers of the atoms to be permuted are shown in red. The dynamic bond "1>0" in CGRs corresponds to a broken single bond, and "0>1" corresponds to a formed single bond. Notice, that the chemical distances (the number of dynamic bonds in CGR) are equal to 2 for both erroneously and correctly mapped reactions.

Figure 3 .

 3 Figure 3. Example of a fixing rule application.Here, an addition-elimination reaction was mapped by ChemAxon toolkit which mixed up the nitrogen atoms N6 and N14. To fix this problem, (a) a CGR is generated, and (b) isomorphic embedding of an erroneous RCext is performed. (c) AAM atomic labels of the nitrogens N6 and N14 in the product are permutated according to the fixing rule. Finally, (d) a new CGR is generated, and the correct RCext is embedded to check if the AAM correction is done successfully. Here, the dynamic bond "1>0" corresponds to a broken single bond, and "0>1" corresponds to a formed single bond.

Figure 4 .

 4 Figure 4. Average CPU time per reaction for different AAM tools .

Figure 5 .

 5 Figure 5. Performance of different AAM tools and consensus strategies.

Figure 6 .

 6 Figure 6. Example of a bromination reaction wrongly mapped by the NameRXN tool. In this reaction, one of the bromine atoms of thionyl bromide (number 19 or 21) replaces the OH group in [1-(4-fluorophenyl)-1,3-dihydro-2-benzofuran-5-yl]methanol. However, in the product, the tool assigns number 23 to the Br atom, which formally corresponds to a fictive reactant missed in the reaction equation.

Figure 7 .

 7 Figure 7. Example of Diels-Alder reaction with erroneous (top) and correct (bottom) atom -to-atom mappings.Here, the reaction itself is on the left side, and the corresponding CGR is on the right side. The symbol "n>m" corresponds to the transformation of bond order n to m[33] . Notice that the chemical distance corresponding to erroneous mapping (CD=5) is lower than that for correctly mapped reaction (CD=6).

Table 1 .

 1 AAM benchmarking results.

	AAM tool or	# correctly	Accuracy,% #reactions	# "unsuitable"	# erroneous	# erroneous
	consensus	mapped		correctly	reactions a	AAM	with	AAM with the
	strategy	reactions		mapped by the		different CD b	same CD c
				given tool only
	ChemAxon	1304	70.45	5	0	339		208
	Indigo	712	38.47	4	518	457		164
	NameRXN	622	33.60	9	873	113		243
	RXNMapper	1550	83.74	49	9	208		84
	RDTool	1411	76.23	0	30	250		160
	New RDTool	1410	76.18	0	9	269		163
	MIN	1316	71.10	0	8	395		132
	MAX	1369	73.96	0	8	311		163
	MIXTURE	1347	72.77	0	8	339		157
	Consensus 1 d	1494	80.71	-e	0	289		68
	Consensus 2	1532	82.77	-	0	246		73
	Consensus 3	1556	84.06	-	0	220		75
	Consensus 4	1554	83.95	-	0	222		75
	Consensus 5	1549	83.68	-	0	221		81
	Consensus 6	1422	76.82	0	8	262		159

a The mapping tool can not map a given reaction b Chemical Distance for the erroneously mapped reaction differs from that with correct AAM.

Full Paper


Tingting Fu

Janina Oetjen

Manuel Chapelle

Alexandre Verdu

Matthias Szesny

Arnaud Chaumot

Davide Degli-Esposti

Olivier Geffard

Yohann Clément

Arnaud Salvador

In situ isobaric lipid mapping by MALDI-ion mobility separation-mass spectrometry imaging

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Conclusions

• MALDI TIMS Imaging leads to richer molecular information content.

• Isobaric and isomeric lipids were mapped to specific locations in a tiny crustacean sample.

• The combination of in-situ MALDI MS/MS and ESI-PASEF is a powerful tool to retrieve lipid identifications with different distributions.

ASMS 2020, TP

Methods

Results

The workflow for MALDI-MSI sample preparation and data analysis is outlined in Figure 1. Fresh frozen female gammarid was sectioned at 12 µm and thaw mounted onto an ITO slide (Sigma-Aldrich). After drying, sections were sprayed with 10 mg/ml DHB in ACN/H 2 O/TFA (70:30:0.1, v/v/v) using a TM sprayer (HTX Technologies, Chapel Hill, NC, USA). Tissues were imaged using the following parameters if not indicated otherwise: m/z range: 100-1000, 400 shots, 10 kHz laser frequency, pitch: 20 µm. For ion mobility separation, ions were separated and eluted in the second part of the dual TIMS device using a ramp time of 300 ms and a 1/K 0 range of 0.6-1. The highly diverse chemical structures of lipids make their analysis directly from biological tissue sections extremely challenging. Here we report the in-situ mapping and identification of lipids in a freshwater crustacean Gammarus fossarum using MALDI mass spectrometry imaging (MSI) in combination with an additional separation dimension using ion mobility spectrometry (IMS) [1]. The high-resolution trapped ion mobility spectrometry (TIMS) allowed efficient separation of isobaric/isomeric lipids showing distinct spatial distributions. The structures of the lipids were further characterized by MS/MS analysis. It is demonstrated that MALDI MSI with mobility separation is a powerful tool for distinguishing and localizing isobaric/isomeric lipids. • Good sensitivity at high lateral resolution.

Gammarus fossarum

References 1. Fu et al., doi: 10.1002/jms.4531 2. https://metlin.scripps.edu 3. http://lipidmaps.org/ 4. Picache et al., Chem. Sci., 2019, 10, 983 Contact: tingting.fu@univ-lyon1.fr; sophie.ayciriex@univ-lyon1.fr

  8. Imaging data processing was performed with SCiLS Lab version 2020a. Ion mobility data was visualized with TIMS data viewer. Assignment of the ions were achieved by interrogating open source databases including METLIN [2], Lipid Maps [3] and CCS Compendium[4].

Fig

  Fig. 3: MALDI-IMS-MSI of female G. fossarum tissue section with 20 µm pixel size.

Fig

  Fig. 1: MALDI-MSI workflow

Fig. 4 :

 4 Fig. 4: High resolution MALDI imaging resolves fine tissue structures in a very small sample. The three-color overlay of ions at m/z 794.606 (red), m/z 772.525 (green) and m/z 375.065 (blue) illustrates the specific distributions of different chemical species in the organs of female G. fossarum.

Table 1 :BFig. 5 :

 15 Fig. 4


Thierry Aubert
email: thierry.aubert@lpmi.uhp-nancy.fr

Omar Elmazria

Badreddine Assouar

Eloi Blampain

Ahmad Hamdan

Damien Genève

Sylvain Weber

Investigations on AlN/Sapphire Piezoelectric Bilayer Structure for High-Temperature SAW Applications

Keywords: high temperature, AlN, sapphire, SAW, langasite

This paper explores the possibility to use AlN/Sapphire piezoelectric bilayer structure for hightemperature surface acoustic wave (SAW) applications. In order to determine the temperature stability of AlN, homemade AlN/Sapphire samples are annealed in air atmosphere for 2 to 20 hours at temperatures going from 700 to 1000°C. Ex situ X-ray diffraction measurements reveal that the microstructure of the thin film is not impacted below 1000°C. Otherwise, ellipsometry and secondary ion mass spectroscopy investigations attest that AlN/Sapphire is reliable up to 700°C. Beyond this temperature, both methods evidence ongoing surface oxidation of AlN. Furthermore, Pt/Ta and Al interdigital transducers are patterned at the surface of the AlN film. The resulting SAW devices are characterized up to 500°C and 300°C respectively, showing reliable frequency response and a large, quasi-constant, temperature sensitivity, with a first order temperature coefficient of frequency around -75 ppm/°C. Besides, between room temperature and 300°C, both electromechanical coupling coefficient K 2 and propagation losses increase, so that the evolution of delay lines insertion losses with temperature strongly depends on the length of the propagation path.

I. INTRODUCTION

Surface acoustic wave (SAW) devices are key components of telecommunication systems for more than thirty years. They also constitute a very promising solution for sensors applications, as they are very sensitive to the external environment, while being passive. This last property is particularly interesting in harsh environments, such as high temperatures up to 1000°C, where this technology offers the possibility to achieve wireless sensors [START_REF] Fachberger | Properties of radio frequency Rayleigh waves on langasite at elevated temperature[END_REF].

The choice of the piezoelectric substrate, which is an essential part of SAW devices, is critical to realize such sensors. Indeed, there are very few piezoelectric materials capable of withstanding high temperatures above 500°C [START_REF] Aubert | Wireless and batteryless surface acoustic wave sensors for high temperature environments[END_REF]. The most promising results have been obtained so far with langasite (La3Ga5SiO14 or

LGS). In particular, Pereira da Cunha et al. have successfully tested this material for 5½ months at 800°C [START_REF] Peireira Da Cunha | High temperature stability of langasite surface acoustic wave devices[END_REF]. Moreover, Hornsteiner et al. have measured a SAW signal on an LGS-based device up to 1085°C [START_REF] Hornsteiner | Surface acoustic wave sensors for high temperature applications[END_REF].

Last but not least, good quality LGS-wafers are commercially available. However, LGS has also some drawbacks. In particular, its acoustic propagation losses increase drastically with frequency and temperature, which restrains its use in ISM wireless mode to the 434 MHz-band, excluding the more promising 2.45 GHz one [START_REF] Fachberger | Properties of radio frequency Rayleigh waves on langasite at elevated temperature[END_REF]. Furthermore, several authors have observed at high temperature a significant increase in the surface roughness of LGS substrates, which onset temperature depends on the gas environment [START_REF] Bardong | Influence of packaging atmosphere on the durability of high-temperature SAW sensors[END_REF]- [START_REF] Schulz | Diffusion-related implications for langasite resonator operation[END_REF]. This phenomenon is related to chemical transformations happening in the vicinity of the surface of LGS crystals. Thus, Bardong et al. observed Ga losses from 650°C under light vacuum environment (pressure = 0.4 mbar), using energy dispersive X-ray spectroscopy [START_REF] Bardong | Influence of packaging atmosphere on the durability of high-temperature SAW sensors[END_REF]. With the help of secondary ion mass spectroscopy (SIMS) measurements, Shulz et al. have quantified those losses in air atmosphere. The quantity of Ga atoms that would leave the crystal through its surface during a period of 100 days at 750°C is in the order of 10 16 per centimeter square [START_REF] Schulz | Diffusion-related implications for langasite resonator operation[END_REF]. These results show that it is mandatory to find new piezoelectric materials which could overcome the limitations of langasite for ultra-high temperature (over 800°C) and ultra-high frequency (over 1 GHz) SAW applications.

Aluminium nitride (AlN) thin films could be one of these materials. It has been identified as a potential candidate for high-temperature SAW applications more than 10 years ago [START_REF] Hornsteiner | Surface acoustic wave sensors for high temperature applications[END_REF], [START_REF] Damjanovic | Materials for high temperature piezoelectric transducers[END_REF]. Indeed, this material is chemically stable up to 1040°C under ultra-high-vacuum conditions [START_REF] Ambacher | Thermal stability and desorption of Group III nitrides prepared by metal organic chemical vapor deposition[END_REF], 1200°C in N2/H2 mixture [START_REF] Lin | Oxidation behavior of AlN films at high temperature under controlled atmosphere[END_REF] and is expected to oxidize in air atmosphere in the 800°C-range [START_REF] Chaudhuri | Thermal oxidation of single crystalline aluminum nitride[END_REF]. Furthermore, showing the highest SAW velocity among the piezoelectric materials, it is particularly suited to carry out devices operating at high frequencies such as 2.45 GHz [START_REF] Caliendo | Theoretical and experimental investigation of gigahertz-band, temperaturecompensated electromechanical coupling configurations based on AlN films[END_REF]. Eventually, some high-quality AlN films show acoustic propagation losses well lower than those of LGS (Tab. I). Despite all these potentialities, AlN has hardly been experimentally studied as piezoelectric material for high-temperature SAW applications. We have recently demonstrated that IDT(Pt/Ta)/AlN/Sapphire SAW structure withstands temperatures as high as 900°C for short-time exposure of 30 minutes [START_REF] Thierry | Surface acoustic wave devices based on AlN/Sapphire structure for high temperature applications[END_REF]. Now, the goal of the present study is to investigate the reliability of the AlN/Sapphire bilayer structure for longer time under such environment. (0001) sapphire wafers [START_REF] Aubert | Highly textured growth of AlN films on sapphire by magnetron sputtering for high temperature surface acoustic wave applications[END_REF]. X-ray diffraction (XRD) measurements as well as transmission electron microscopy revealed that the (002) AlN films are stress-free and highly-textured out of the plane, with a full width at half-maximum (FWHM) of the (002) rocking-curve (RC) as low as 0.26°. Both methods also evidenced that these films are textured in the plane as well. The value of the grain size, deducted from XRD measurements and atomic force microscopy, is close to 40 nm. The latter method also gave access to the surface roughness of the films which value, equal to 8 Å (root mean square), is suitable for SAW applications.

To study the temperature stability of the AlN/Sapphire structure in air atmosphere, the 3" wafers were cut in small pieces (20 x 20 mm). The resulting samples were annealed at temperatures between 700 and 1000°C, for periods going from 2 to 20 hours in a tubular furnace with open ends. The effects of high-temperature exposure were studied by doing XRD (θ-2θ and rocking-curve scan modes, PANalytical X'Pert Pro MRD), ellipsometry (Horiba Jobin Yvon Uvisel) and SIMS (primary ions: Cs + , Cameca IMS 7F) measurements.

Furthermore, to perform interdigital transducers (IDTs) withstanding high-temperature conditions, 10 nm thick tantalum as adhesion layer and 100 nm thick platinum films as electrodes, were deposited by sputtering method on some AlN/Sapphire substrates [START_REF] Aubert | Behavior of platinum/tantalum as interdigital transducers for SAW devices in high-temperature environments[END_REF]. SAW delay lines were then fabricated by photolithography and ion beam etching. The devices were constituted by two identical IDTs, with 50 pairs of fingers and a wavelength λ fixed to 16 μm. The IDT center spacing, acoustic aperture, and finger width to space ratio were 175λ, 70λ, and 1:2, respectively. The SAW propagation path was along the Sapphire Xdirection. These devices were electrically characterized up to 500°C in air atmosphere using a Network Analyzer (Agilent PNA 5230A, Santa Clara, CA) and RF prober station (Signatone S-1160, Gilroy, CA) equipped with Signatone thermal probing system S-1060 Series so that it can withstand temperatures up to 600°C. However, we restricted our characterization to 500 °C due to the limitation imposed by the use of cooled RF probes (Z-Probe Suss Microtech). Finally, to get the evolution of the electromechanical properties of the AlN/Sapphire bilayer structure with temperature, conventional 100 nm thick Al IDTs were used to minimize the influence of the electrodes on the devices response. These measurements were conducted up to 300°C, considering the electromigration phenomena that aluminium undergoes from this temperature [START_REF] Thiele | Platinum and palladium high-temperature transducers on langasite[END_REF]. Namely, the propagation losses behavior was studied from delay lines with the same design as the one described above, except that the IDT center spacing was comprised between 55λ and 110λ in that case. Some synchronous resonators were also fabricated to access the electromechanical coupling coefficient K 2 of the structure, via the following equation [START_REF] Hines | A simple transducer equivalent circuit parameter extraction technique[END_REF]:

where Go is the conductance value at the central frequency f0, C the capacitance between IDT finger pair and N the number of finger pairs in the IDT. To do so, the wavelength was reduced to 6 μm so that the SAW is more confined in the piezoelectric film, leading to a higher value of K 2 , and by the way to a better accuracy in the determination of its evolution.

III. RESULTS AND DISCUSSION

A. Temperature stability of AlN in air atmosphere

On the figures 1 and 2, the RC of the XRD (002) AlN peak and the c-lattice parameter of the AlN films are respectively plotted as function of the annealing duration for temperatures going from 700 to 1000°C.

The c-lattice parameter was deduced via XRD measurements, from the Bragg angle θ of the (002) AlN peak, using the Bragg formula [START_REF] Tiggelaar | Stability of thin platinum films implemented in high-temperature microdevices[END_REF]. Its evolution is related to the one of the stress in the film. One can see on these two figures that neither the texture nor the stress of the AlN layer is modified by high-temperature exposure, provided the temperature is lower than 1000°C. These results evidence the physical stability of the AlN/Sapphire bilayer structure under such environments. On the other side, the intensity of the XRD (002) AlN peak decreases all along the annealing process. Figure 3 illustrates this phenomenon in the case of a 2h-long annealing at 1000°C. As this cannot be related to a decline of the film crystalline quality (texture or grain size), its origin has to be found in the shrinking of the AlN amount, related to its progressive oxidation from the surface [START_REF] Chaudhuri | Thermal oxidation of single crystalline aluminum nitride[END_REF]. However, XRD does not allow to evidence the formation of an aluminium oxide overlayer below 1000°C, likely because this layer is amorphous if the annealing temperature is too low [START_REF] Yue | SIMS study on the initial oxidation process of AlN ceramic substrate in the air[END_REF]. Figure 4 summarizes the evolution of the intensity of the (002) AlN peak with time for temperatures going from 700 to 1000°C.

As expected, the kinetics of this oxidation reaction strongly depends on the annealing temperature. It follows an exponential law which leads to an asymptotic regime in a matter of hours. This kind of behavior is related to the fact that the diffusion of the oxygen molecules through the oxide overlayer, responsible for further oxidation, is getting more and more difficult as the thickness of the latter increases [START_REF] Kang | Periodic oxide breakdown during oxidation of AlN/Sapphire(0001) films[END_REF]. After 20h of annealing, the asymptotic regime is reached except for the highest temperature in this test, i.e. 1000°C.

One can observe that after such duration, the intensity of the (002) AlN peak is less than 10% of its initial value at 1000°C, around 50% at 900°C and more than 90% at 800 and 700°C. As this intensity is directly related to the AlN amount [START_REF] Kingetsu | In situ diffractometry of cristobalite formation during hightemperature oxidation of SiC films[END_REF], it means that the order of magnitude of the oxide overlayer thickness is respectively of 1 μm, 500 nm and some tens of nanometers. Thus, for the lower temperatures, the variations in thickness of the oxide overlayer are likely too small to be correctly quantified by ex situ XRD measurements, as attests the important dispersion of the results at 700°C visible on figure 4. of the oxide layer thickness from the one of the FWHM of this SIMS peak [START_REF] Benninghoven | Secondary Ion Mass Spectrometry[END_REF]. One can observe the presence of a thin native oxide layer at the surface of the as-deposited AlN films. Its thickness slightly increases after 2 hours of annealing at 700 and 800°C but remains close to that of the native one, which is not the case at higher temperatures, thereby confirming XRD measurements. The most important result appears after 20 hours of annealing. Whereas the thickness of the oxide layer continues to slightly increase at 800°C, it seems that the initial minor oxidation passivates AlN films from further aggression at 700°C. A more quantitative estimation of the oxide layer thickness was obtained by performing ellipsometric measurements. To do so, one has to make experimental optical parameters and calculated ones from dedicated softwares (in our case, DeltaPsi2 from Horiba Jobin Yvon) fit together. The accuracy of the method depends mainly on the correct modeling of the sample structure. In our case, reliable results with high correlation coefficients were obtained by simply describing annealed samples as composed by 3 layers one above the other: amorphous Al2O3/AlN/sapphire. The oxide layer thickness deducted from this method for the same samples as the ones investigated by SIMS can be seen on figure 7. The correlation between the results given by both methods is very good. In particular, the main information provided by SIMS method is fully confirmed by ellipsometry: the 10 nm thick native oxide layer is growing a little bit by 7-8 nm during the first 2h of annealing at 700°C, and then is stable during longer annealing process at this temperature. As this stability is not observed at higher temperatures, 700°C seems the upper limit regarding the stability of AlN/Sapphire at high temperature in air atmosphere.

B. SAW measurements

The experimental frequency-temperature law of the Pt/Ta/AlN/Sapphire SAW delay lines is plotted on figure 8. The investigated structure is characterized by a large and quasi-constant sensitivity to temperature, with first-order temperature coefficient of frequency (TCF) values going from -68 ppm/°C at 20°C to -80 ppm/°C at 500°C, making it well suited for temperature sensor applications. Note that as the relative thickness h/λ of the AlN layer (where h is the thickness of the AlN film) is quite limited here, with a value smaller than 0.1, the SAW is mainly located in the sapphire substrate. Thus, this thermal behavior can be mainly attributed to the latter. Figure 9 shows the time-gated frequency response (S21 magnitude) of the same delay line at different temperatures between the ambient and 500°C. The main observation is the quasi-constant increase in the insertion losses at a rate close to 1 dB/100°C. Regarding delay lines, such an improvement should result in a decrease of insertion losses with temperature. However, as it can be seen on figure 11, the evolution of the latter strongly depends on the length of the propagation path. In fact, whereas a significant improvement (insertion losses diminution) can be seen for the shortest delay line between room temperature and 300°C, one can observe a relative stability for intermediate ones and a slight deterioration in the case of the longer one. This means that the propagation losses also increase when the temperature is rising and that this phenomenon can counteract the expansion of K 2 if the delay line is long enough. In the case of the Pt/Ta/AlN/Sapphire delay line investigated up to 500°C, the propagation length was even more important with a value of 125λ, which explains this increase of the insertion losses at a non-negligible rate of 1 dB/100°C. The results obtained in the first part of this study suggest that the AlN/Sapphire bilayer structure can be used as piezoelectric material for SAW sensors operating in air atmosphere up to 700°C, as the oxidation process seems quite limited, which is not the case at higher temperatures. This statement will have to be confirmed by in situ SAW measurements conducted up to this temperature, in order to determine the impact of this very thin additional oxide overlayer. Anyway, better stability of AlN is expected in environments less aggressive than air, such as those containing small amounts of oxygen, where temperatures as high as 1000°C could be likely reached [START_REF] Lin | Oxidation behavior of AlN films at high temperature under controlled atmosphere[END_REF]. Moreover, the AlN/Sapphire bilayer structure could have an enormous hidden potential for high-temperature SAW applications, as the IDTs could be placed below the AlN film. Thus, AlN/IDT/AlN/Sapphire or AlN/IDT/Sapphire configurations are both promising solution leading to protect both IDTs from agglomeration phenomena and active piezoelectric area of AlN from oxidation, or other chemical aggression, not to mention that this would significantly improve the electromechanical coupling coefficient K 2 of the structure [START_REF] Caliendo | Theoretical and experimental investigation of gigahertz-band, temperaturecompensated electromechanical coupling configurations based on AlN films[END_REF]. The results presented in the second part of this paper, dedicated to SAW measurements on devices working around as it shows a large, quasi-constant, temperature sensitivity up to at least 500°C. We also observed that both K 2 and propagation losses are increasing with temperature so that the evolution of the insertion losses of delay lines based on this structure depends mainly on their length. Finally, though this study was conducted at relatively low frequency (around 350 MHz), these results are also encouraging for highfrequency applications (i.e. mainly in the ISM band at 2.45 GHz), for which AlN/Sapphire structure shows great interest over LGS. In that case, if the AlN thickness was kept around 1.3 μm, its relative value h/λ would reach 0.6, leading to a wave more confined into the piezoelectric film and thus to an electromechanical coupling coefficient K 2 likely not too far from the maximum value for this structure, i.e. 0.3%. Concerning the TCF, such an increase in the relative thickness of AlN should result in a reduction of its value to 20-30 ppm/°C at room temperature [START_REF] Tsubouchi | Zero-Temperature-Coefficient SAW devices on AlN epitaxial films[END_REF], which is still suitable for temperature sensor applications. On the other side, it is not possible to predict if in that case, the frequency-temperature law will still show a quasi-linear behavior in the 20-500°C range as the thermal behavior of pure AlN is not known. At last, a larger operating frequency means as usual more propagation losses. But, as Uehara et al.

showed, if the AlN film is of high quality with very low roughness (some Å) and low dislocation density, then the AlN/Sapphire structure could show well-contained propagation losses, even lower than those of lithium niobate [START_REF] Uehara | Fabrication of 5-GHz-Band SAW filter with atomically-flat-surface AlN on Sapphire[END_REF]. The only remaining question will be then the quantitative evolution of these propagation losses with temperature.

Fig. 1 .

 1 Fig. 1. FWHM of the rocking-curve of AlN (002) as a function of the annealing duration

Fig. 3 .Fig. 4 .

 34 Fig. 3. θ-2θ XRD scan mode of as-prepared and 2h annealed AlN/Sapphire samples

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. SIMS depth profile signal of an AlN/Sapphire sample after a 2h long annealing process. A thin conductive Pt layer was sputtered at the surface to reduce charging effects.

Fig. 8 .Fig. 9 .Fig. 10 .

 8910 Fig. 8. Temperature dependence of the synchronous frequency of the delay line

Fig. 11 .

 11 Fig. 11. Delay line length influence on the insertion losses evolution with temperature

  TAB. I. COMPARISON BETWEEN SAW PROPERTIES OFLGS AND AlN/SAPPHIRE BILAYER STRUCTURE[START_REF] Aubert | Wireless and batteryless surface acoustic wave sensors for high temperature environments[END_REF]

	Material	SAW velocity (m/s)	Electromechanic al coupling coefficient K 2 (%)	Acoustic propagation losses at 1 GHz (mdB/λ)
	Langasite (0°, 138.5°, 26.6°)	2700	0.4	4
	(002)AlN/(0001)Sapphire (khAlN = 2)	5600	0.3	0.7
		II. EXPERIMENTAL
	1.3 μm thick (002) AlN thin films were epitaxially grown by reactive magnetron sputtering onto 3"


Vitaly Solov'ev
email: solovev-vp@mail.ru

Alexandre Varnek
email: varnek@unistra.fr

A N Frumkin

Thermodynamic radii of lanthanide ions derived from metal-ligand complexes stability constants

Keywords: Linear free energy relationships, stability constant prediction, complexation thermodynamics, complexation similarity, complexes of metal ions with organic ligands, lanthanide metal ion, organic ligand

Here, we describe two sets of thermodynamic ionic radii (r1 and r2) of lanthanide ions derived from the analysis of large set of the stability constants log1 and log2 of the ML and ML2 complexes of lanthanide ions M (from Ce 3+ to Lu 3+ ) with organic ligands (L) in water. It has been demonstrated that the stability constants of two metals Mi and Mj with a given ligand L are related by simple equations log1 j = r1 i /r1 j. log1 i and log2 j = r2 i /r2 j. log2 i which formally correspond to purely electrostatic interactions between spherical cations with organic molecule. Predictive performance of these equations was assessed in 5-fold cross-validation procedure. The standard deviation (s) of predictions varies from 0.3 to 1.0 for log1 and from 0.4 to 1.2 for log as a function of the difference of Shannon ionic radii (rij) of Mi and Mj: s1 = 0.27 + 4.25rij for log1 and s2 = 0.39 + 5.16rij for log. The new radii r1 and r2 steadily decrease across the Ln series, i.e., they follow the same trend as Shannon effective radii of lanthanide ions. The calculations were performed using experimental data on 2854 log1 values for 445 organic ligands and 947 log2 values for 156 organic ligands in the complexes with 13 metals.

Introduction

Lanthanides (Ln) are widely used in modern innovative technologies including lanthanidedoped nanomaterials, superconductors, and magnetic resonance imaging [START_REF] Speghini | Synthesis of Upconverting Nanomaterials: Designing the Composition and Nanostructure[END_REF][START_REF] Bunzli | Luminescence Bioimaging with Lanthanide Complexes[END_REF][START_REF] Chen | Lanthanide-Doped Luminescent Nanomaterials. From Fundamentals to Bioapplications[END_REF][START_REF] Cotton | A Perspective on Lanthanide Chemistry[END_REF][START_REF] Li | Luminescent Rare Earth Complexes as Chemosensors and Bioimaging Probes[END_REF] . Industrial applications of rare earths have developed in metallurgy, magnets, ceramics, electronics, chemical, optical, medical, nuclear technologies. The published data on world rare earth production demonstrates that annual production exceeded 120000 metric tons of rare earth oxides [START_REF] Gupta | Extractive Metallurgy of Rare Earths[END_REF] . The chemical properties of lanthanides are very similar, because when moving from La to Lu, the structure of the external electron shell of the Ln 3+ cations remains and electrons fill only the inner 4f layer [START_REF] De Bettencourt-Dias | Introduction to Lanthanide Ion Luminescence[END_REF] . Similar chemical properties lead to the big difficulty in lanthanide separation processes. Currently, pure compounds of these metals are produced by their extraction separation in two-phase system from aqueous solution into organic solvent using organic molecules (ligands, extractants), capable to form with Ln 3+ ions complexes soluble in organic solvents [START_REF] Krishnamurthy | Extractive Metallurgy of Rare Earths[END_REF][START_REF] Alyapyshev | Recovery of Minor Actinides from High-level Vastes: Modern Trends[END_REF][START_REF] Ustynyuk | N'-Dialkyl-N,N'-diaryl-1,10-phenanthroline-2,9-dicarboxamides as donor ligands for separation of rare earth elements with a high and unusual selectivity. DFT computational and experimental studies[END_REF] . Extraction technologies of the same type are used for the separation of lanthanides and actinides in the nuclear fuel cycle [START_REF] Alyapyshev | Recovery of Minor Actinides from High-level Vastes: Modern Trends[END_REF] . The development of highly selective ligands is of great importance because differences in the stability constants of the complexes have been used for effective separation of lanthanides [START_REF] Gupta | Extractive Metallurgy of Rare Earths[END_REF] .

For this purpose, quantum-chemical modeling is widely used and makes it possible to significantly reduce the number of compounds to be synthesized and studied experimentally [START_REF] Ustynyuk | Quantum chemical modelling of extraction separation of minor actinides and lanthanides: the state of the art[END_REF][START_REF] Lan | Recent advances in computational modeling and simulations on the An(III)/Ln(III) separation process[END_REF] .

Computational costs are significantly lower for a combination of quantum chemistry and machine learning methods [START_REF] Montavon | Machine learning of molecular electronic properties in chemical compound space[END_REF][START_REF] Pereira | Machine Learning Methods to Predict Density Functional Theory B3LYP Energies of HOMO and LUMO Orbitals[END_REF][START_REF] Solov'ev | Predictive Models for HOMO and LUMO Energies of N-Donor Heterocycles as Ligands for Lanthanides Separation[END_REF] or machine learning methods using experimental data regarding the stability constants of the complexes and various quantitative parameters of extraction. Different machine learning methods (MLR, SVM, ANN, kNN and others) and various types of descriptors were involved for regression modeling of stability constants of the Ln 3+ -ligand complexes in solution [START_REF] Solov'ev | Classification of Metal Binders by Naive Bayes Classifier on the Base of Molecular Fragment Descriptors and Ensemble Modeling[END_REF][START_REF] Solov'ev | New Approach for Accurate QSPR Modeling of Metal Complexation: Application to Stability Constants of Complexes of Lanthanide Ions Ln3+, Ag+, Zn2+, Cd2+ and Hg2+ with Organic Ligands in Water[END_REF][START_REF] Solov'ev | Complexation of Mn2+, Fe2+, Y3+, La3+, Pb2+, and UO22+ with Organic Ligands: QSPR Ensemble Modeling of Stability Constants[END_REF][START_REF] Tetko | Benchmarking of Linear and Non-Linear Approaches for Quantitative Structure-Property Relationship Studies of Metal Complexation with Organic Ligands[END_REF][START_REF] Dioury | QSPR Prediction of the Stability Constants of Gadolinium(III) Complexes for Magnetic Resonance Imaging[END_REF][START_REF] Svetlitski | QSPR Modelling of Lanthanide-Organic Complex Stability Constants[END_REF][START_REF] Kiani-Anbouhi | Prediction of the complexation stabilities of La3+ ion with ionophores applied in lanthanoid sensors[END_REF][START_REF] Kiani-Anbouhi | Application of QSPR for prediction of the complexation stabilities of Sm(III) with ionophores applied in lanthanoid sensors[END_REF][START_REF] Qi | Correlation Analysis of the Structures and Stability Constants of Gadolinium(III) Complexes[END_REF] .

Most of these models were developed on small training sets and validated on only one specifically selected test set [START_REF] Solov'ev | Classification of Metal Binders by Naive Bayes Classifier on the Base of Molecular Fragment Descriptors and Ensemble Modeling[END_REF] ; only some modeling workflows involved a cross-validation step and ligands of various organic classes [START_REF] Solov'ev | New Approach for Accurate QSPR Modeling of Metal Complexation: Application to Stability Constants of Complexes of Lanthanide Ions Ln3+, Ag+, Zn2+, Cd2+ and Hg2+ with Organic Ligands in Water[END_REF][START_REF] Solov'ev | Complexation of Mn2+, Fe2+, Y3+, La3+, Pb2+, and UO22+ with Organic Ligands: QSPR Ensemble Modeling of Stability Constants[END_REF][START_REF] Tetko | Benchmarking of Linear and Non-Linear Approaches for Quantitative Structure-Property Relationship Studies of Metal Complexation with Organic Ligands[END_REF][START_REF] Solov'ev | The complexation of metal ions with various organic ligands in water: prediction of stability constants by QSPR ensemble modelling[END_REF][START_REF] Solov'ev | QSPR Ensemble Modelling of the[END_REF][START_REF] Solov'ev | QSPR ensemble modelling of alkaline-earth metal complexation[END_REF][START_REF] Solov'ev | Stability Constants of Complexes of Zn2+, Cd2+, and Hg2+ with Organic Ligands: QSPR Consensus Modeling and Design of New Metal Binders[END_REF] .

Linear free energy relationships (LFERs) are used as a simple method for stability constant predictions. The basis of LFER is the existence of a linear relationship between the logarithm of the equilibrium constants (free energies) for two reactions at the varied reagent or parameter of equilibria [START_REF] Dimmock | Approaches to predicting stability constants[END_REF][START_REF] Hancock | Approaches to Predicting Stability Constants. A Critical Review[END_REF] , e.g. for the series of ligands L, the presence of a regression equation between the stability constants (logK) of the complexes of two metal ions Mi and Mj: logKMiL = a + b . logKMjL. Choppin established LFER of the stability constants of aminocarboxilate ligands between Am 3+ and Nd 3+ 31 .

The LFER of the stability constants of the trivalent actinides (Pu, Am, Cm, Bk and Cf) and the similarly sized lanthanides were used to estimate the stability constants for 91 amino(poly)carboxilates [START_REF] Uhnak | Prediction of trivalent actinide amino(poly)carboxylate complex stability constants using linear free energy relationships with the lanthanide series[END_REF] . Simple regression equations without the intercept term a (logKMiL = b . logKMjL) are discovered between Co 2+ and Ni 2+ , Co 2+ and Cu 2+ , Ni 2+ and Cu 2+ for various 272, 221 and 339 organic ligands, respectively [START_REF] Solov'ev | QSPR Ensemble Modelling of the[END_REF] , and between Al 3+ and Fe 3+ , In 3+ and Fe 3+ , Th 4+ and VO 2+ , for small series of the ligands [START_REF] Solov'ev | The complexation of metal ions with various organic ligands in water: prediction of stability constants by QSPR ensemble modelling[END_REF] .

In this paper, we analyze the relationships log i = a + b . log j between the stability constants (log) of the complexes of ith and jth lanthanide ions M 3+ (from Ce 3+ to Lu 3+ ) with diverse organic ligands (L). The calculations were performed on the stability constants log1  logK1 and log2 respectively for the complexation reactions M 3+ + L = (M 3+ )L and M 3+ + 2L = (M 3+ )L2 in water at 298 K and an ionic strength 0.1 M. All pair combinations of lanthanide ions were studied excepting Pm 3+ . Using developed LFER for every pair of lanthanides, two sets of "thermodynamic" ionic radii of lanthanide ions r1 for (M 3+ )L and r2 for (M 3+ )L2 were derived for an accurate prediction of stability constants for ith lanthanide ion using stability constants for jth lanthanide ion. That let us to predict the stability constants log1 and log2 of numerous potential complexes for different lanthanide/ligand combinations. The standard deviations of the predictions are defined by a difference of ionic radii of two metal ions involved in the calculations.

Methods

Data Sets

The datasets used in this study comprise the logarithms of the stability constants (log1  logK1 and log2) for the ML and ML2 complexes of 13 lanthanide ions (Ce 3+ , Pr 3+ , Nd 3+ , Sm 3+ , Eu 3+ , Gd 3+ , Tb 3+ , Dy 3+ , Ho 3+ , Er 3+ , Tm 3+ , Yb 3+ , Lu 3+ ) with various organic ligands in water. The experimental log1 and log2 values were critically selected from the IUPAC Stability Constants Database (SC DB, version 5.86, Academic Software) [START_REF] Pettit | IUPAC Stability Constants Database[END_REF] at standard temperature 298 K and an ionic strength I = 0.1 M. A few log1 and log2 values were corrected to specified temperature and ionic strength using the van't Hoff equation and the Davies equation, respectively, by the procedures included in SC DB. If several values of log1 (log2) were available for a particular ligand, the most recent data or the data consistent with respect to different experimental methods were chosen.

The number of ligands in each metal dataset is given in Table 1. The smallest subset used in this study corresponds to Lu 3+ ligands (53 log2 values), while the largest one to Gd 3+ ligands (308 log1 values). The ranges of stability constant values vary from 0.2 (Nd 3+ ) to 26.4 (Yb 3+ ) for log1 and from 2.8 (Er 3+ ) to 33.1 (Eu 3+ ) for log2 (Fig. 1). The datasets contain 445 (156) organic ligands and 2854 (947) log1 (log2) values. The ligands include derivatives of cyclic and acyclic polydentate molecules of various organic classes. As a rule, the organic ligand has several functional groups. Groups of the ligands can include amino, mercapto and hydroxy derivatives of carboxylic acids; various aminocarboxylates; different aminoacids and their oligomers; derivates of phosphoric and phosphinic acids; alkyl-and aminophosphonic acids; cyclic and acyclic polydentate ligands with the terminal carboxy groups separated by various cyclic or acyclic spacers; various sulfonic acids; ternary amines with phosphono and carboxy groups; mono-and dipodands of ternary amines; phenols, nitroso and amino derivatives of phenols; adenosine and guanosine derivatives; crown-ethers, thia-, and aza-crowns with neutral and acidic lariat groups, cryptands; pyridines; phenanthroline and 1,3-dione derivatives; cyclodextrins, etc.

(Supplementary Material: Table SM1 and the Structure Data Files: Ln_logBeta1_445L.SDF and Ln_LogBeta2_156L.SDF).

The external test set includes 12 polydentate cyclic and acyclic organic ligands and 31 log1 values. The datasets comprise recent experimental log1 for lanthanides Ce 3+ , Eu 3+ , Gd 3+ , Ho 3+ , Lu 3+ and Sm 3+ (Table SM10).

Table 1. The complexation of lanthanide ions (M 3+ ) with diverse organic ligands (L) in water at 298 K and an ionic strength 0.1 M. The numbers of studied ligands for the equilibria M 3+ + L = (M 3+ )L (log1) and M 3+ + 2L = (M 3+ )L2 (log2

Computational Procedures

2D structures of the ligands, names of lanthanide ions and corresponding log1 and log2 values resulted from searching in SC DB were converted into the Structure Data Files (SDF) which served as the initial data for the LFER modeling. The SDF manager EdiSDF [START_REF] Solov'ev | Structure-Property Modelling of Complex Formation of Strontium with Organic Ligands in Water[END_REF][START_REF] Varnek | Substructural Fragments: An Universal Language to Encode Reactions, Molecular and Supramolecular Structures[END_REF][START_REF] Solov'ev | Editor of Structure -Data Files[END_REF][START_REF] Varnek | ISIDA -Platform for Virtual Screening Based on Fragment and Pharmacophoric Descriptors[END_REF] was used to prepare the datasets and to build LFERs. An EdiSDF option was used for calculations and graphical analysis of linear regressions between each pair of quantitative data fields of SDF.

The program LanComSim (Lanthanide Complexation Similarity) [START_REF] Solov'ev | Lanthanide Complexation Similarity[END_REF] was created for predicting the stability constants log1 and log2 on the basis of thermodynamic ionic radii of lanthanides (see sections 3.2 and 3.4).

The external 5-fold cross validation (5CV) [START_REF] Solov'ev | QSPR Models on Fragment Descriptors In Tutorials in Chemoinformatics[END_REF] and external test set were applied to validate relationships with using thermodynamic ionic radii for predicting the stability constants. Predictive performance of the relationships has been estimated by root mean squared error RMSE and squared determination coefficient 𝑅 𝑑𝑒𝑡 2 estimated on the combination of external test sets in 5CV or on the external test set: Here 𝑦 𝑖 and ŷ 𝑖 are, respectively, experimental and predicted values of log1 (log2) for ith ML (ML2) complex, n is the number of data points (the stability constant pairs, ligands), < 𝑦 > is the mean of experimental log1 (log2) value. The Pearson's correlation coefficient (R), the standard deviation (s) and the Fisher criterion (F) are used to describe regression equations.

𝑅𝑀𝑆𝐸 = [ 1 𝑛 ∑(𝑦 𝑖 -ŷ 𝑖 ) 2 𝑛 𝑖=1 ] 1/2 𝑅 𝑑𝑒𝑡 2 = 1 -∑(𝑦 𝑖 -ŷ 𝑖 ) 2 / ∑(𝑦 𝑖 -< 𝑦 >) 2 𝑛 𝑖=1 𝑛 𝑖=1 log1

Results and Discussion

From Linear Free Energy Relationships (LFER) to Thermodynamic Ionic Radii

The LFER method was applied to predict the stability constants of the complexes of lanthanide ions (M 3+ ) with diverse organic ligands (L). The modeling was performed on the stability constants log1 and log2 of the complexations M 3+ + L = (M 3+ )L and M 3+ + 2L = (M 3+ )L2 in water at 298 K and an ionic strength 0.1 M, respectively. For the series of the diverse sets of the ligands L (Table 1, Supporting information: Table SM1), the regression equations 1 and 2 between the stability constants of the complexes were built for all pairs of lanthanide ions Mi and Mj.

log1 j = a + b . log1 i (1) log2 j = a + b . log2 i (2)
Here a and b are coefficietnts of the equation for the ith and jth lantanide ions and series of the ligands. For 13 lanthanide ions from Ce 3+ to Lu 3+ excepting Pm 3+ , 156 regression equations (eqs) 1 (2) were built for all pairs of metal ions (Tables SM2 andSM3). Eqs have the reasonable squared correlation coefficient (R 2 ) and standard deviation (s): 0.960 ≤ R 2 ≤ 0.999 and 0.16 ≤ s ≤ 0.97 for the 156 eqs. 1, and 0.890 ≤ R 2 ≤ 0.999 and 0.17 ≤ s ≤ 1.9 for the 156 eqs. 2 (Tables SM2 andSM3).

For eqs, the number of the data points (n) changes from 111 (Ce 3+ vs. Tm 3+ ) to 223 (Nd 3+ vs. Sm 3+ ) for log1 and from 29 (Ce 3+ vs. Lu + ) to 73 (Sm 3+ vs. Gd 3+ ) for log2.

The calculated intercept term a and its standard deviation sa vary respectively from -0.24 to 0.37 and from 0.02 to 0.18 in eq. 1, from -0.41 to 1.2 and from 0.27 to 0.73 in eq. 2. The a and sa values are very similar (Tables SM2 andSM3) and compatible with random errors of log1 (log2) [START_REF] Pettit | IUPAC Stability Constants Database[END_REF] , therefore we built simple linear regressions 3 and 4 without the intercept term a:

log1 j = b . log1 i (3) log2 j = b . log2 i (4)
The standart deviations of eqs. 1 and 3 can be compared by F-test = s 2 eq 1/s 2 eq 3 [START_REF] Muller | Tafeln der mathematischen Statistik[END_REF] . F-test varies from 0.90 to 1.01, and the critical value of F155, 155; 0.95 40 is 1.32 at the 0.10 significance level, we can conclude that seq 1 = seq 3. Therefore, the simple eq. 3 can be used instead of the eq. 1 predicting log1 (Table SM4). The similar conclusion is true concerning eq. 4 and the log2 predictions (Table SM5). The simple eqs. 3 and 4 have as good statistical parameters as the eqs. 1 and 2: 0.960 ≤ R 2 ≤ 0.999 and 0.16 ≤ s ≤ 0.97 for the 156 eqs. 3, 0.879 ≤ R 2 ≤ 0.999 and 0.17 ≤ s ≤ 2.0 for the 156 eqs. 4 (Tables SM4 andSM5).

The simple equation 3 (4) means, that the coefficients b of the ensemble of the 156 eqs.

(Tables SM4 andSM5) are interdependent, where one is the product of two others, and this indicates high complexation thermodynamics similarity of the lanthanide ions. Every lanthanide ion can be described by one characteristic r1 (r2) for the ML (ML2) complex: log1 j /log1 i = r1 i /r1 j = b (log2 j /log2 i = r2 i /r2 j = b) according to eq. 3 (4). The system of the 156 equations r1 i /r1 j = b (r2 i /r2 j = b) for log1 (log2) was solved to calculate 13 unknown characteristics r1 (r2) applying nonlinear regression analysis using the simplex algorithm [START_REF] Himmelblau | Applied Nonlinear Programming[END_REF] . Two sets of "thermodynamic" ionic radii r1 and r2 for the lanthanide ions were derived (Table 2) for an accurate prediction of the stability constant log1 j (log2 j ) for jth lanthanide ion and the given ligand using the stability constant for ith lanthanide ion and eq. 5 (eq. 6). R 2 and s for eqs. 3 (4) and 5 ( 6) are very similar:

for most metal pairs, the difference is in the fourth (R 2 ) and third (s) decimal places (Tables SM4 -SM7

).

log1 j = r1 i /r1 j log1 i (5)

log2 j = r2 i /r2 j. log2 i (6)
The radii r1 and r2 vary from 0.983 to 1.087 (r1) and from 0.970 to 1.160 (r2). From Ce 3+

to Lu 3+ , r1 decreases by 9.5%. However, this decrease is not monotonous. The radius r1 is reduced by 7.1% from Ce 3+ to Sm 3+ , becomes about constant (Sm 3+ -Gd 3+ ) and then declines slightly by 1.7% (Tb 3+ -Lu 3+ ). The radius r2 decreases by 16.3% (Ce 3+ -Lu 3+ ), the increment is also not monotonous with a knee at Sm 3+ (Fig. 2). The new radii r1 and r2 steadily decrease across the Ln series, i.e., they follow the same trend as Shannon effective radii 42 of lanthanide ions. The radius r2 corresponds to the Shannon effective ionic radius rCN8 for ion coordination number (CN) 8 (Fig.

2):

r2 = (-0.027  0.047) + (1.028  0.045)rCN8, n = 13, R 2 = 0.979, s = 0.008

From Ce 3+ to Eu 3+ , the radius r1 corresponds to the Shannon ionic radius rCN7 (Table SM11), then it decreases to a lesser extent across the Ln series as compared with rCN7, and from Er 3+ to Lu 3+ it is very close to the Shannon radius rCN8 (Fig. 2). The thermodynamic radius r2 of the ML2 complexes containing two ligands naturally correcponds to the Shannon radius with a higher coordination number compared to the thermodynamic radius r1 of the ML complexes. Ionic radii of Ln 3+ cations decrease steadily across the Ln series, hence heavier Lns often have lower coordination numbers than lighter Lns [START_REF] Mills | Ligand Design in Modern Lanthanide Chemistry[END_REF] , but heavier Lns have higher stability of the complexes than lighter Lns according to the thermodynamic radii r1 and r2. A similar relationship was observed in the series of transition metals [START_REF] Tsebrikova | Complexation of the new tetrakis[methyl(diphenylphosphorylated)] cyclen derivative with transition metals: First examples of octacoordinate zinc(II) and cobalt(II) complexes with cyclen molecules[END_REF] : coordination number of metal ion in complexes decreases in the Co 2+ L -Ni 2+ L -Cu 2+ L order, but log1 increases for the same ligand.

Predicting the Stability Constants by Thermodynamic Ionic Radii

The fivefold external cross validation (5CV) [START_REF] Solov'ev | QSPR Models on Fragment Descriptors In Tutorials in Chemoinformatics[END_REF] was used to evaluate the predictive performance of the relationships 5 and 6. In this procedure, the entire dataset of the ligands was split into five non-overlapping pairs of training and test sets. Each training set covers 4/5th of the initial dataset while the related test set covers the remaining 1/5th. The "thermodynamic" radius r1 (r2) fitted on the i-th training set were applied to the corresponding test set, thus, all the ligands of the initial data set are predicted including various combinations of lanthanide ion pairs. The relationship 5 (6) performs well in cross-validation comparing 25474 (8002) experimental and predicted log1 (log2) pairs (Fig. 3). Squared determination coefficient 𝑅 𝑑𝑒𝑡 2 is high: 0.987 (r1) and 0.982 (r2). Root mean squared error RMSE is 0.57 (r1) and 0.76 (r2) which is compatible with the systematic errors in known experimental log1 (log2) values [START_REF] Pettit | IUPAC Stability Constants Database[END_REF] .

An accuracy of the predictions depends on similarity of metal ion pairs in eqs. 5 and 6. The standard deviations of the linear regressions 5 and 6 on the base of the radii r1 and r2 (Tables SM6 andSM7) increase as a function of the difference of Shannon ionic radii (rij) of Mi and Mj (Fig. 4), and the predicted log1 and log2 by eqs. 5 and 6 are described with well defined standard deviations s1 and s2, respectively: Here r is the Shannon effective ionic radius for coordination number 8 [START_REF] Shannon | Revised Effective Ionic Radii and Systematic Studies of Interatomie Distances in Halides and Chaleogenides[END_REF] . The regression equations 7 and 8 demonstrate that the estimated standard deviation s increases with the rij difference from 0.3 to 1.0 for log1 and from 0.4 to 1.2 for log2. The thermodynamic radii r1 and r2 give a possibility to avoid labour-consuming experimental measurements of the stability constants for every lanthanide ion. For instance, it is enough to evaluate the stability constants of the complexes of metal ions such as Nd 3+ , Tb 3+ and Tm 3+ . This ensures low standard deviations of the predicted constants for the 10 remaining lanthanides. For 10 of 445 ligands (Tables SM8 andSM9 in SM), the pairs of the experimental log1 (log2) values for two different metal ions are incompatible with the regression equation 5 (6). The difference between the experimental and predicted stability constants varies from 2.6 to 5.8 (log1) and from 2.9 to 8.4 (log2) which significantly exceeds the standard deviation estimated by eqs. 7 and 8. This could be explained by either exceptional selectivity of these ligands to some specific lanthanide ions or, most likely, inaccurate experimental estimations of stability constants.

s1 = 0.27 + 4.25rij (
It is known, that the logarithm of a ratio of the stability constants for two cations log(1 j /1 i ) is used as a complexation selectivity measure [START_REF] Solov'ev | Complexation of Phosphoryl-Containing Mono-, Bi-and Tri-Podands with Alkali Cations in Acetonitrile. Structure of the Complexes and Binding Selectivity[END_REF] . The equation 5 (6) demonstrates that high complexation selectivity log(1 j /1 i ) = (1 -r1 j /r1 i )log1 i (log(2 j /2 i ) = (1 -r2 j /r2 i )log2 i ) among the lanthanide ions can be reached by ligands, the complexes of which are characterized very high stability constants. In the case of extraction Ln separation in two-phase system from aqueous solution into organic solvent, the high stability constant is important in the organic phase [START_REF] Alyapyshev | Recovery of Minor Actinides from High-level Vastes: Modern Trends[END_REF][START_REF] Ustynyuk | Quantum chemical modelling of extraction separation of minor actinides and lanthanides: the state of the art[END_REF][START_REF] Ivanova | Effect of the 2-[(diphenylphosphoryl)metoxyphenyl]diphenylphosphine oxide derivative structures on the extraction and ion-selective properties toward rare-earth elements[END_REF] and easier to implement, because the stability of the complexes in many organic solvents is much higher than in water [START_REF] Pettit | IUPAC Stability Constants Database[END_REF][START_REF] Solov'ev | Supramolecular chemistry .60. Solvent effects on crown ether complexations[END_REF] . In addition, the formation of the ML2 complexes is preferable, as a rule, due to the higher solubility in organic solvent compared to the ML complexes.

External Test Set

The relationship 5 was also tested on the log1 predictions for the external test set of the complexes of recently synthesized organic ligands. The set includes 12 ligands [START_REF] Kálmán | Synthesis, Potentiometric, Kinetic, and NMR Studies of 1,4,7,10-Tetraazacyclododecane-1,7-bis(acetic acid)-4,10-bis(methylenephosphonic acid) (DO2A2P) and its Complexes with Ca(II), Cu(II), Zn(II) and Lanthanide(III) Ions[END_REF][START_REF] Försterová | Thermodynamic study of lanthanide(III) complexes with bifunctional monophosphinic acid analogues of H4dota and comparative kinetic study of yttrium(III) complexes[END_REF][START_REF] Baranyai | The effects of intramolecular H-bond formation on the stability constant and water exchange rate of the Gd(III)-diethylenetriamine-N'-(3-amino-1, 1-propylenephosphonic)-N,N,N'',N''-tetraacetate complex[END_REF][START_REF] Aguilar-Pérez | Phosphate Ester Hydrolysis by Hydroxo Complexes of Trivalent Lanthanides Stabilized by 4-Imidazolecarboxylate[END_REF][START_REF] Marques | 13-and 14-membered macrocyclic ligands containing methylcarboxylate or methylphosphonate pendant arms: Chemical and biological evaluation of their 153Sm and 166Ho complexes as potential agents for therapy or bone pain palliation[END_REF] for which the experimental data were not presented in SC DB. The ligands are 12-, 13-and 14-membered azacrowns containing acetate or/and methylenephosphonate pendant groups (Table SM10 in SM). The log1 values were predicted for the 31 complexes at 298 K and an ionic strength 0.1 M excepting 4-imidazolecarboxylic acid for which an ionic strength was 0.05 M. The predicted log1, pred demonstrate reasonable agreement (Rdet 2 = 0.989 and RMSE = 0.74) with the observed log1, exp (Fig. 5).

Predictor of stability constants

A predictor of the stability constants log1 (log2) of the ML (ML2) complexes of lanthanide ions M 3+ with diverse organic ligands (L) in water was created on the base of the relationships 5 -8 and the "thermodynamic" ionic radii r1 and r2 (Table 2). The predictor is publically available on the server http://vpsolovev.ru/programs/ [START_REF] Solov'ev | Lanthanide Complexation Similarity[END_REF]

Conclusions

The LFER method was used to study relationships between the stability constants log1 (log2) of the ML (ML2) complexes for different pairs of lanthanide ions M 3+ with the same ligands (L). It was shown that the LFERs log1 j = a + b . log1 i (log2 j = a + b . log2 i ) derived from the analysis of large set of the stability constants are able to predict the stability constants with high accuracy compatible with the systematic errors in known experimental log1 (log2) values. Here i and j are indexes of pairs of lanthanide ions from Ce 3+ to Lu 3+ for the complexes M i L and M j L (log1), M i L2 and M j L2 (log2) in water at 298 K and an ionic strength 0.1 M. Then we demonstrated that the intercept term a is not significant, and simple ratios can be used, from which it follows that every lanthanide can be described by the thermodynamic ionic radius r1 (r2) for the ML (ML2) complexes, and the stability constants of two metal ions Mi and Mj with a given ligand L are related by the simple equations log1 j = r1 i /r1 j. log1 i and log2 j = r2 i /r2

= a + b . log1 i (log2 j = a + b . log2 i )
for the stability constant log1 (log2) of the complexation M 3+ + L = (M 3+ )L (M 3+ + 2L = (M 3+ )L2) of lanthanide ions M 3+ i and M 3+ j with the organic ligands (L) in water at 298 K and an ionic strength 0.1 M. Table SM4 (SM5) contains the statistical parameters of 156 regression equations log1 j = b . log1 i (log2 j = b . log2 i ) for the stability constant log1 (log2) of the complexation M 3+ + L = (M 3+ )L (M 3+ + 2L = (M 3+ )L2) of lanthanide ions M 3+ i and M 3+ j with the organic ligands. Table SM6 (SM7) contains the statistical parameters of 156 regression equations log1 j = r1 i /r1 j. log1 i (log2 j = r2 i /r2 j. log2 i ) for the stability constant log1 (log2) of the complexation M 3+ + L = (M 3+ )L (M 3+ + 2L = (M 3+ )L2) of lanthanide ions M 3+ i and M 3+ j with the organic ligands. Table SM8 (SM9) includes the several pairs of the experimental log1 (log2) values which are incompatible with the equations log1 j = r1 i /r1 j. log1 i (log2 j = r2 i /r2 j. log2 i ). Table SM10 comprises experimental and predicted stability constants log1 for the external test set. Table SM11 contains effective ionic radii of the lanthanide ions for different coordination numbers.

The Structure Data Files Ln_logBeta1_445L.SDF and Ln_LogBeta2_156L.SDF comprise experimental data: 2D structures of ligands, 2854 log1 values for 445 organic ligands and 947 log2 values for 156 organic ligands for the ML and ML2 complexes of 13 lanthanide ions (Ce 3+ , Pr 3+ , Nd 3+ , Sm 3+ , Eu 3+ , Gd 3+ , Tb 3+ , Dy 3+ , Ho 3+ , Er 3+ , Tm 3+ , Yb 3+ , Lu 3+ ) in water at standard temperature 298 K and an ionic strength I = 0.1 M. Fig. 4 The standard deviations (s1 and s2) of the linear regressions log1 j = r1 i /r1 j. log1 i (a) and log2 j = r2 i /r2 j. log2 i (b) versus the difference of Shannon ionic radii (rij) of lanthanide ion pairs M 3+ i and M 3+ j.

FIGURE CAPTIONS

Fig. 5 Experimental versus predicted stability constants log1 for the external test set.

TABLES CAPTIONS Table 1. The complexation of lanthanide ions (M 3+ ) with diverse organic ligands (L) in water at 298 K and an ionic strength 0.1 M. The numbers of studied ligands for the equilibria M 3+ + L = (M 3+ )L (log1) and M 3+ + 2L = (M 3+ )L2 (log2).

Table 2. Thermodynamic ionic radii r1 and r2 of lanthanide ions Ln 3+ .

Fig. 1

 1 Fig. 1 Ranges of experimental log1 (a) and log2 (b) for the ML (a) and ML2 (b) complexes. The left and right boundaries of the boxes indicate the 25 th and 75 th percentiles; a line within the box marks the median. Error bars of the boxes indicate the 90 th and 10 th percentiles.

Fig. 2 .

 2 Fig. 2. Thermodynamic ionic radii r1 (a) and r2 (c), Shannon effective ionic radii (Å) for coordination numbers 8 (b), 7 (d) and 6 (e) of lanthanide ions Ln 3+ .

Fig. 3

 3 Fig. 3 Predictive performance of the relationships 5 and 6 in 5CV procedure: predicted versus experimental stability constant log1 (a) and log2 (b).

Fig. 4

 4 Fig.4 The standard deviations (s1 and s2) of the linear regressions log1 j = r1 i /r1 j. log1 i (a) and log2 j = r2 i /r2 j. log2 i (b) versus the difference of Shannon ionic radii (rij) of lanthanide ion pairs M 3+ i

Fig. 5

 5 Fig. 5 Experimental versus predicted stability constants log1 for the external test set.

Fig. 1

 1 Fig. 1 Ranges of experimental log1 (a) and log2 (b) for the ML (a) and ML2 (b) complexes. The left and right boundaries of the boxes indicate the 25 th and 75 th percentiles; a line within the box marks the median. Error bars of the boxes indicate the 90 th and 10 th percentiles.

Fig. 2 .Fig. 3

 23 Fig. 2. Thermodynamic ionic radii r1 (a) and r2 (c), Shannon effective ionic radii (Å) for coordination numbers 8 (b), 7 (d) and 6 (e) of lanthanide ions Ln 3+ .

  ).

	0	5	10	15	20	25	5	10	15	20	25	30
	no. metal		the number of ligands
		ion		log1		log2
	1	Ce 3+		198		66
	2	Pr 3+		240		83
	3	Nd 3+		277		88
	4	Sm 3+		249		83
	5	Eu 3+		239		79
	6	Gd 3+		308		95
	7	Tb 3+		186		82
	8	Dy 3+		221		77
	9	Ho 3+		202		64
	10 Er 3+		190		61
	11 Tm 3+		157		58
	12 Yb 3+		198		58
	13 Lu 3+		189		53

Table 2 .

 2 Thermodynamic ionic radii r1 and r2 of lanthanide ions Ln 3+ .

	no. metal ion	r1	radius	r2
	1	Ce 3+	1.087		1.160
	2	Pr 3+	1.064		1.136
	3	Nd 3+	1.047		1.113
	4	Sm 3+	1.015		1.066
	5	Eu 3+	1.012		1.061
	6	Gd 3+	1.015		1.051
	7	Tb 3+	1.000		1.038
	8	Dy 3+	0.998		1.029
	9	Ho 3+	0.997		1.018
	10	Er 3+	0.995		1.013
	11	Tm 3+	0.989		1.000
	12	Yb 3+	0.984		0.997
	13	Lu 3+	0.983		0.970

  Predictive performance of these equations was assessed in 5-fold cross-validation procedure. The standard deviation (s) of predictions varies from 0.3 to 1.0 for log1 and from 0.4 to 1.2 for log as a function of the difference of Shannon ionic radii (rij) of Mi and Mj: s1 = 0.27 + 4.25rij for log1 and s2 = 0.39 + 5.16rij for log. The new radii r1 and r2 steadily decrease across the Ln series, i.e., they follow the same trend as Shannon effective radii of lanthanide ions. The calculations were performed using experimental data on 2854 log1 values for 445 organic ligands and 947 log2 values for 156 organic ligands in the complexes with 13 metals. Using the developed program LanComSim we predicted unknown 2931 (1081) of log1 (log2) values of the potential complexes.Table SM1 contains the names of the studied organic ligands. Table SM2 (SM3) contains the statistical parameters of 156 regression equations log1 j

	correspond to purely electrostatic interactions between spherical cations with organic molecule.
	log1, pred
	log1, pred = 0.13 + 0.995 log1, exp
	n = 31
	Rdet 2 = 0.989
	RMSE = 0.74
	Supplementary Material and Software
	log1, exp
	j. log2 i which formally

A predictor 38 of the stability constants log1 (log2) of the ML (ML2) complexes of lanthanide ions M 3+ with diverse organic ligands (L) in water is freely available as the program LanComSim (Lanthanide Complexation Similarity) under Windows operating systems.

Conflict of Interest

The authors declare that they have no conflict of interest.