Yuta Miyazaki

Dr Christophe Kahlfuss

Dr Ayumu Ogawa

Dr Takashi Matsumoto

Dr Jennifer A Wytko

Prof. T Koji Oohora

Takashi Hayashi
email: thayashi@chem.eng.osaka-u.ac.jp

Dr Jean Weiss
email: jweiss@unistra.fr

Prof Jean-Pierre Sauvage

Hayashi

CuAAC in a Distal Pocket: Metal Active-Template Synthesis of Strapped-Porphyrin [2]Rotaxanes

Keywords: active template, cycloaddition, CuAAC, rotaxane, porphyrinoid

In modern supramolecular chemistry, the copper catalysed alkyne-azide cycloaddition (CuAAC), [1] better known as the click reaction, has provided access to a large variety of elaborate, but otherwise difficult to prepare molecular structures. [2][3][4][5] Several detailed analyses of the role of copper in the mechanism of the cycloaddition are available in the literature. [7][8][9] It is now commonly admitted that two metal cations are involved in the cycloaddition step. [10] This cooperative process is supported by experimental evidence in recent reports. [3,9,11] The exact role played by the second metal is not clearly defined and can range from activation, due to Lewis acidic character, to coordination of the zwitterionic form of the azide reagent. The cooperative and simultaneous activation of both the azide and the alkyne by two distinct metal centres is compatible with both a catalytic [6] or stoichiometric use of copper in click reactions. The latter use strongly relates to the concept of an active metal template in the click synthesis of rotaxanes that has been developed [12] and employed by several groups. [12][13][14][15] In this approach, the active participation of two metal centres in the mechanism of the click reaction was highlighted by the use of copper-copper [12] or copper-nickel [13] in the formation of rotaxanes. This report demonstrates that the use of a ditopic ligand able to bind both an active copper template and a coordinatively unsaturated metal leads to both active and organising template effects and provides access to mechanically linked porphyrin architectures. Supramolecular architectures built on a phenanthroline (phen)-strapped porphyrin 1(H2), which is available on gram scale, have been the focus of our group for almost two decades. This ditopic ligand has led to several species in which the phen binds and stabilizes a copper(I) ion, [16,17] whereas the porphyrin binds a variety of metals, thus offering the possibility of coordinating an exogenic substrate on either the open face, opposite to the strap, or within the phen pocket. [18,19] As confirmed by both NMR and Xray studies, the ligand is bound within the phen pocket in inclusion complexes of 1(Zn) with imidazole. [16] Thus, the possible occurrence of a click reaction within the phen strap of 1(Zn) was anticipated due to the isomorphic character of 1,2,3-triazole and imidazole, and an assistance by the coordination of an azido reagent to the zinc of the porphyrin was also expected. Scheme 1. Synthesis of 4(H2). Reactions were performed in the presence of 0.96 equivalents of Cu I , and subsequently treated with EDTA and TFA.

Treatment of the zinc porphyrin 1(Zn) with a slightly substoichiometric amount of Cu(MeCN)4PF6 yielded the zinc-copper complex in situ. As shown in Scheme 1, reaction in basic medium of the alkyne 2 [12a] and azide 3, [20] bearing bulky stoppers commonly used in mechanically interlocked structures, [21] afforded mixtures of rotaxane structures still containing either copper, zinc or both metals. Literature reports of the active template method for rotaxane synthesis have established the possible extrusion of copper resulting from sterically demanding scaffoldings. [12][13][14][15] To avoid tedious purification of mixtures containing partially metalated species, treatment with ethylenediaminetetraacetate (EDTA) to fully remove the copper was performed. Because this treatment initiated a partial demetalation of the zinc, the latter was then completely removed upon trifluoroacetic acid (TFA) treatment and only the free base rotaxane species 4(H2) was isolated in various yields (Entries 1-7) and characterised. The corresponding nickel and zinc rotaxanes 4(Ni) and 4(Zn) were later obtained by metalation of 4(H2) with nickel and zinc salts, respectively (see SI).

To demonstrate the essential role played by the metal in the porphyrin core, the nature of the metal in the porphyrin was investigated (Table 1). Due to the previously reported fast disproportionation of a copper(I) phen-strapped porphyrin and its rapid conversion to the copper(II) porphyrin, synthetic attempts involving the free base porphyrin 1(H2) and the copper(I) phen strap complex were precluded. [17] Even though zinc(II) salts have been reported to catalyse azide-alkyne cycloadditions, [22] the use of 1(Zn) alone in the absence of Cu(MeCN)4PF6 did not produce rotaxane or the free dumbbell under the same reaction conditions, emphasizing the expected role of copper(I) as the active metal template.

The use of 1(Cu), which weakly binds an axial ligand, resulted in a significantly lower yield of rotaxane (Entry 8) compared to that obtained with 1(Zn). This result shows that the coordination on the metalloporphyrin is of utmost importance for a productive dual template effect. Indeed, the 16% yield obtained in the presence of 1(Cu) is representative of the oxophilic nature of hard porphyrin-Cu(II) complexes, [START_REF] Cremers | For the quantification of Cu-pyridine interaction (G~0) see[END_REF] whereas the soft porphyrin-Zn(II) complex has a clearly enhanced affinity for nitrogen-containing axial ligands and affords up to 50% yield of rotaxane (Entry 5). Finally, when the nickel porphyrin 1(Ni), which is inert towards the axial binding of exogenic ligands, was employed, no rotaxane was detected (Entries 9-11). Therefore, the presence of a second, coordinatively unsaturated metal in the porphyrin core is necessary for the success of the click reaction.

Table 1. Yields of rotaxane 4(H2) obtained by click-reaction of 2 and 3 in the presence of 1(M) and 0.96 eq. of Cu(MeCN)4PF6 Entry M Base [a] Yield (%) [b] 1 Zn none n.d. The coordination of the azide on the zinc in the porphyrin core of 1(Zn) is in agreement with the observed influence of the base on the efficacy of the CuAAC reaction (Entries 1-7). Thus, the nature of the base used in the click reaction was also studied. The yield of the rotaxane increased with increasing basicity, from picoline to 2,6-lutidine and DIPEA, emphasizing the role of the base in the deprotonation of the terminal alkyne and stressing the necessary weak binding of a nitrogen ligand on the open face of the porphyrin (Entries 1-7). The geometry of the azide group leaves a choice for coordination of the zwitterionic intermediate on the open face (Figure 1 A) or on the strapped-side (Figure 1 C or D). However in the presence of the phen-Cu-acetylide, bases, such as picoline and DIPEA, have no other choice than to bind to To rule out possible geometric effects influencing the efficacy of the cycloaddition, the dimensions of the strapped porphyrin scaffolds of [1(Zn-Cu)], [18] 1(Ni), 4(Ni) and 4(H2) were compared in their respective crystal structures (see Figure 2 andS12). The dimensions of the hollow pyramidal space between the porphyrin plane and the phen strap are nearly identical in all four structures; therefore, geometric considerations cannot explain the discrimination observed in the reactivity of 1(Zn) and 1(Ni) in the presence of copper(I) for the templated 1,3-cycloaddition of the azide and alkyne stoppers. Thus, only the nature of the metal in the porphyrin and its ability to assist the templated process are responsible for the formation of the rotaxane. In this regard, the classical Sauvage macrocycle, derived from 2,9diphenylphenanthroline pentaethylene-glycol, did not produce any rotaxane [12b] by the active metal template method, confirming the necessary assistance by the metal in the porphyrin core observed in our work. Thus, two distinct and complementary roles are played by each metal. Copper is the active metal template, as previously described in works by Leigh [12,13] and Goldup. [14] Zinc plays the role of an organising template according to Sauvage's definition [24] and is responsible for a productive positioning of the azide in the cavity defined by the phen strap, probably as shown in Figure 1 D. At this stage, the necessary competing binding of a labile, encumbered, axial ligand on the porphyrin's open face, which favours the productive binding of the azide in the strap, limits the efficiency of the zinc organising template effect to 50%. Current work is in progress to explore variation of the base or the ditopic nature of the macrocycle to enhance the productive binding of the azide in the strap.

Single crystals of 4(H2) and 4(Ni) suitable for X-ray diffraction analysis were obtained by vapour diffusion of acetonitrile into a saturated solution of 4(H2) in dichloroethane/acetonitrile (1/1) and vapour diffusion of acetonitrile a saturated CH2Cl2 solution of 4(Ni) (see SI for details). In the solid state structure of 4(H2) (Figure 2a) the propyl chain's carbon atoms and the porphyrin plane are separated by distances ranging from 3.3 to 3.6 Å, which are close to van der Waals distances. Weaker C-H/ interactions probably contribute to the additional stabilisation of the propylene chain threaded through the phen strap. The triazole develops - stacking interactions with the porphyrin ring at a distance of ca. 3.7 Å and the dumbbell's phenyl group on the left-hand side in Figure 2a is stacked above the edge of the porphyrin ring at a distance of ca. 3.3 Å. The latter stabilising interaction was previously observed for an inclusion complex of 2phenylimidazole in 1(Zn). [25] In the crystal structure of 4(Ni) in Figure 2b, the weak oxophilic interaction, denoted by the 3.0 Å distance between the nickel core and the phenol ether, is assisted by two additional noncovalent interactions. The stacking of the left-hand phenyl group and the porphyrin edge occurs again, with distances ranging from 3.6 to 3.9 Å, in spite of the slight distortion of the porphyrin macrocycle. Furthermore, the protons located on the first and third carbons of the propylene chain are at rather short distances (2.3 and 2.7 Å) from the phen's nitrogen atoms.

These structural features are persistent in solution, as confirmed by 2D 1 H NMR experiments (see supporting information). In all three rotaxanes, displacements of chemical shifts observed for different segments of the dumbbell's signals indicated the positioning of each dumbbell segment in the phenstrapped porphyrin (Figure S7). In 4(Zn), the triazole was clearly located in the phenanthroline-porphyrin pocket, as shown by the deshielding of the triazole proton that is hydrogen bonded to the phenanthroline. [14e] The binding of the triazole to zinc was also detected in the UV-visible spectrum of 4(Zn) where the Soret and Q bands appeared at 428 and 559 nm rather than at 418 and 546 nm in the spectrum of 1(Zn) (see SI). The red-shifted absorptions of 4(Zn) are consistent with reported observations of triazole coordinated to a zinc porphyrin. [26] In conclusion, this work confirms the role played by a second metal ion, zinc(II), in the 1,3-dipolar cycloaddition mechanism and opens the way to a "tandem active metal template" approach to rotaxanes that cannot be obtained by a "single active metal template" method. Although this example uses a ditopic porphyrin ligand, the approach could be extended to a wide range of ditopic ligands. In the particular case of easily accessible phenanthrolinestrapped porphyrins, this strategy is currently being explored for the straightforward preparation of mechanically linked multiporphyrin scaffolds such as photochemical dyads and triads.

Experimental Section

Procedures, full characterisation, 1 H NMR NOESY spectra, crystal structure, data and tables for 1(Ni), 4(Ni) and 4(H2) are available in Supporting Information.

Figure 1 .

 1 Figure 1.Proposed tandem actions of a Zn organising template and a Cu active metal template: The competitive binding of the azide on the metallic porphyrin core and the axial base is ruled by their respective affinities for the central metal.

  zinc on the porphyrin's open face, as schematically shown in Figure1B. An assistance of the click reaction by coordination of the azide to a second metal requires binding of the azide on the strappedside of the porphyrin, as proposed in Figure1C or D. Azide binding can be oriented to the strapped-side if a fifth ligand, in fast exchange, is already present on the open face of the zinc porphyrin. Picoline strongly binds to zinc but only on the open face of the porphyrin due to steric requirements; therefore, this base is in slow exchange and clearly competes with the binding of azide on the strapped-side (Entry 2). With DIPEA, which is a stronger yet non-nucleophilic base that weakly binds to zinc, the temporary residence of this base on the open face of the zinc porphyrin favours the insertion of the azide on the strapped-side (Entry 5). With 2,6-lutidine, the steric hindrance in the vicinity of the nitrogen atom prevents binding on either side of the porphyrin, thus allowing both the efficient coordination of the azide on the strapped-side and the non-productive azide binding on the open face (Entry 3). To try to take advantage of the reactivity of a strong base and the weak binding of a bulky base on the open face, stoichiometric amounts of both diazabicycloundecene (DBU) and DIPEA were employed; however, the 45% yield was slightly lower than when using only DIPEA (Entry 7)

Figure 2 .

 2 Figure 2. Solid state structures of rotaxanes 4(H2) (a) and 4(Ni) (b). Solvents were removed for clarity. Significant distances are discussed in the text.

Acknowledgements

We thank Prof. Kazushi Mashima and Dr. Haruki Nagae, Department of Materials Engineering Science, Osaka University, for their help in determining the structure of 1(Ni). The authors are grateful for financial support from the Université de Strasbourg, the CNRS, Osaka University, JSPS (AO), Incoming International Master program of IdEx of the Université de Strasbourg (YM), a JST-ANR grant AI\R.1 4.JTIC.OOO2.O 1 (MECANO) (CK) and the JSPS KAKENHI grant JP15H05804.

Entry for the Table of Contents

COMMUNICATION

CuAAC in a ditopic strapped porphyrin: Two metals cooperate for the click-reaction leading to new [2]rotaxanes derived from a ditopic phenanthroline-strapped porphyrin. In addition to the copper active-template effect, coordination of the azide to the porphyrin's metal core is necessary to generate sterically constrained [2]rotaxanes.


Michael Coeurdassier

Nicolas Bassin

Thibaut Powolny

Christophe Morin

Eve Afonso

Christophe Morin

Damien Bertin

Violaine Champion

Natacha Leclerc

Geneviève Magnon

Dominique Michelat

Jacques Morel

François Rey-Demaneuf

François Richard

Michel Sauret

Bruno Tissot

Prevalence of hematozoan parasites in Red Kite nestlings

Keywords: Haemosporida, Haemoproteus, Leucocytozoon, Plasmodium, raptor Aknowledgements

Blood pathogens may be harmful for birds and thus, represent a potential threat for some populations. To date, only few data are available on the prevalence of hematozoans in the Red Kite, a raptor endemic to the Western Palearctic listed as Near Threatened within its global range. Here, we report the first screening of haemosporidian parasites in Red Kite nestlings from French breeding areas. The blood of 47 nestlings was sampled in three areas in Eastern France and the prevalence of three genera of blood pathogens, Leucocytozoon, Haemoproteus and Plasmodium, was measured based on PCR-assays. Thirteen nestlings were positive to Leucocytozoon sp. (n=13) and Plasmodium relictum (n=1), which corresponds to an overall prevalence of 28%. This prevalence is higher than those previously reported for the Red Kite in other countries. This is the first record of P. relictum, a potential pathogen for birds, in the Red Kite. No relationship was found between the body condition of nestlings and the occurrence of haemosporidians. Our findings suggest that hematozoans could be a limited threat for nestlings but prevalence in adults, a sensitive stage for the population dynamic of Red Kite, remains poorly known and should be further investigated at least in its main breeding areas.

Declarations

Funding: This study was financially supported by the French Ministère de l'Ecologie

Introduction

The Red Kite Milvus milvus is an Accipitridae endemic to the Western Palearctic. It breeds mainly in Germany, France and Spain, which together hold around 75% of the global population (Knott et al. 2010). In these countries, the populations have been estimated to be declining, which led the International Union for the Conservation of Nature to assess the Red Kite as Near threatened (BirdLife International 2020) within its global range while it is currently listed as vulnerable in France (UICN France et al. 2016). The main threats identified are poisoning by anticholinesterasic pesticides and rodenticides, decrease of livestock, homogenization of landscapes and potentially wind turbines (BirdLife International 2020). In both Europe and France, the decline reported during the period 1990-2010 has led to the implementation of a Conservation Plan in which the health status of the Red Kite was surveyed, notably on nestlings that were caught at nests for ringing and wing-marking. Among the health indicators monitored and considering that ecotoxicological issues were identified, chemical concentrations in blood were measured and we showed that 50% of nestlings were exposed to organic compounds and/or metals [START_REF] Coeurdassier | Exposure of nestlings Red Kite to rodenticides[END_REF][START_REF] Powolny | Corticosterone mediates telomere length in raptor chicks exposed to chemical mixture[END_REF]. If infection by pathogens has not been pointed out as a threat for the Red Kite, only few studies have focused on the prevalence of blood parasites and no data are available on the prevalence of hematozoan pathogens in Red Kite nestlings from France. Haemosporidians belonging to three genera of protozoan (Haemoproteus, Leucocytozoon, and Plasmodium) are intracellular blood parasites transmitted to hosts by bloodsucking dipteran vectors, they have been commonly detected in birds including raptors [START_REF] Remple | Intracellular Hematozoa of raptors: a review and update[END_REF]. Haemoproteus and Leucocytozoon are often considered to have little or no pathogenic effects on birds but their presence in some raptor nestlings has been related to subclinical effects and even host death.

Plasmodium, the causative agent of malaria in animals, can be a serious avian pathogen, Plasmodium relictum being the species most commonly associated with virulent disease [START_REF] Remple | Intracellular Hematozoa of raptors: a review and update[END_REF]. In raptors, haemosporidians may also negatively affect health such as body condition and life traits [START_REF] Korpimaki | Blood parasites and reproductive success of Tengmalm's Owls: detrimental effects on females but not on males?[END_REF][START_REF] Hatchwell | The haematozoan parasites of common blackbirds Turdus merula: associations with host condition[END_REF][START_REF] Sol | Parasite mediated mortality and host immune response explain age-related differences in blood parasitism in birds[END_REF][START_REF] Garvin | Blood parasites of nearctic-neotropical migrant passerine birds during spring trans-gulf migration: impact on host body condition[END_REF]. Here, we aim to measure the prevalence of common blood parasites in birds of the genera Haemoproteus, Leucocytozoon and Plasmodium in nestlings of three areas monitored in Eastern France. The relationship between prevalence and body condition of the nestlings was also checked to document whether such infection could affect individuals.

Methods

Study areas and sample collection

Fieldwork was carried out in 2015 from 4 th to 25 th June in Franche-Comté region, Eastern France. Three breeding populations of Red Kite, from three areas (called "Area" factor hereafter) were monitored: Sundgau Belfortain (47°380170'N, 6°510460'E), Drugeon-Remoray basin (46°540220'N, 6°210170'E) and Besançon plateau (47°140350'N, 6°010190'E). Forty-seven chicks (23 males and 24 females) corresponding to 24 nests (called "Nest" factor hereafter) were retained for the study. The sampling was implemented 15 days prior to the expected fledging date, when the chicks were approximately 25-40 days old. Nestlings were removed from the nests by an experienced arborist climber and carried to the ground for data collection. First, the body mass (g) and the wing length (i.e., length from the extremity of the longest primary feather to the wrist of the right wing in mm) were measured. Blood (2 mL) was sampled by trained operators from the brachial vein using lithium-heparin coated syringes. As recommended [START_REF] Fair | Guidelines to the use of wild birds in research[END_REF], the volume of blood collected was <0.5% of the chick's body mass. Blood samples were kept in the dark at 4°C and then stored at -20°C until processed in the laboratory.

The sex of the chicks was determined according to a PCR method [START_REF] Griffiths | A DNA test to sex most birds[END_REF] as described in [START_REF] Powolny | Corticosterone mediates telomere length in raptor chicks exposed to chemical mixture[END_REF]. All the manipulations were performed with governmental authorizations no. 2015093-0007 and 2015093-0023 in the framework of the Red Kite national conservation plan.

Body condition estimate

In the present study, we define body condition as a proxy of energy reserves, mainly fat mass that represents a broad indicator of the health status of nestlings.

According to [START_REF] Labocha | Morphometric indices of body condition in birds: a review[END_REF], body mass of birds is « as good or nearly as good an indicator of fat content, as any of the other condition indices ». Because body mass can be easily and reliably measured, we used it as an estimate of body condition. As we studied growing individuals, difference in body mass could be due to slight differences in specimen age rather than in fat reserve. [START_REF] Mougeot | Breeding biology of the Red Kite Milvus milvus in Corsica[END_REF] showed that the wing length was closely correlated to the age of Red Kite nestlings. Moreover, the relationship they established is almost linear within the range of ages of the kite chicks we studied, thus, we used the wing length as an estimate of chick age. Body condition was calculated as the residuals of the [START_REF] Zuur | Mixed effects models and extensions in ecology with R[END_REF]. The null model was also tested and the most parsimonious model with the lowest Akaike information criterion (AICc) was selected, a difference of AICc > 2 was retained to evidence difference between 2 models [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF]. Normality of model residuals and variance homoscedasticity were verified. All statistical analyses were performed using R 3.4.2. with the packages lme4, pgirmess and lmerTest.

Results

Among the 47 chicks studied, 30 came from the Drugeon-Remoray basin (15 females and 15 males), 14 from the Besançon plateau (eight females and six males) and three from the Sundgau belfortain (one female and two males). Haemosporidian parasites were found in 13 chicks (i.e., 28%), Leucocytozoon was detected in eight males and five females from the three studied areas and Plasmodium in a single male also infected by Leucocytozoon of the Sundgau belfortain population (Fig. 1).

No difference of Leucocytozoon prevalence was evidenced between sex (binomial GLMER, p = 0.33) or between brood size (binomial GLMER, p = 0.66). We found two Leucocytozoon haplotypes, which perfectly matched (i.e. with a 100% identity on 479bp) with L-MILVUS1 (Genbank accession number JN164716) and L-MILANS4 (JN164713) haplotypes described in [START_REF] Pérez-Rodríguez | Molecular characterization of haemosporidian parasites from kites of the genus Milvus (Aves: Accipitridae)[END_REF]. The DNA sequence corresponding to L-MILVUS1 also closely matched with Leucocytozoon toddi (GenBank accession number DQ177272) with a 99% identity.

The DNA sequence amplified in the only sample positive for Plasmodium matched with Plasmodium relictum (GenBank accession number LN835311) with a 100% identity over 462bp. Although Leucocytozoon toddi is a typical parasite of birds [START_REF] Valkiūnas | Identification of Leucocytozoon toddi group (Haemosporida: Leucocytozoidae), with remarks on the species taxonomy of leucocytozoids[END_REF], we cannot exclude that we detected an unreferenced species genetically close, we thus retained Leucocytozoon ssp. hereafter, while we considered that DNA from Plasmodium relictum was detected in the corresponding positive nestling. The only PCR-positive sample for Plasmodium relictum was also positive for Leucocytozoon sp. (L-MILVUS1 haplotype) and we recorded three other individuals with a co-occurrence of both Leucocytozoon haplotypes L-MILVUS1 and L-MILANS4, among which two nestlings from the same brood. This suggests that co-infections can occur relatively frequently, both for interspecific (1/47 individual) and inter-haplotypes (3/47) co-occurrences.

At a same wing length, females are heavier than males (best LMER selected with 'Sex' as factor, p<0.0001, Table 1. Estimate for males = -84.8, df = 31.1, t = -5.45).

Thus, the body condition of the nestlings is presented for both sex separately in Fig. 1. No difference of body condition was detected between infected and non-infected kites whatever the sex (p=0.15, Table 1).

Discussion

Our findings showed that nestlings of Red Kite are infected by haemosporidians, notably Leucocytozoon ssp. and P. relictum, in three French breeding areas. The prevalence of 28% we measured was higher than those commonly reported for Leucocytozoon in Red Kite chicks from German, Spanish, UK or Danish populations while Plasmodium was less frequently found (Table 2). Overall, if we aggregate all the studies reported in Table 2 except our and without distinction of age, the prevalences measured in the Red Kite (n=366) were 3.8%, 0.3% and 0.0% for Leucocytozoon, Haemoproteus and Plasmodium, respectively. Co-infection by two haplotypes of Leucocytozoon corresponding to some of those previously described in Spanish populations of Red Kite (Pérez-Rodriguez et al. 2013) were found in three nestlings of the French populations studied. We also point out that the detection of P. relictum in one single nestling of a French population is the first record for the Red Kite while several haplotypes of Plasmodium were detected in the Black Kite Milvus migrans [START_REF] Pérez-Rodríguez | Molecular characterization of haemosporidian parasites from kites of the genus Milvus (Aves: Accipitridae)[END_REF][START_REF] Nourani | Dynamics of prevalence and distribution pattern of avian Plasmodium species and its vectors in diverse zoogeographical areas -A review[END_REF].

More generally, P. relictum has alrealdy been reported in species belonging to Accipitridae and Falconidae families [START_REF] Bennett | Avian haematozoa: mortality and pathogenicity[END_REF]. It is considered as a serious pathogen for birds and has been associated to virulent diseases and even death for translocated Gyrfalcon Falco rusticolus [START_REF] Remple | Intracellular Hematozoa of raptors: a review and update[END_REF]). However, to our knowledge, no information is available on its harmfulness for Accipitridae and the individual infected by P. relictum we found did not show obvious health troubles.

Prevalence of blood parasites in diurnal raptors is highly variable and may reach 50% or more notably in adults (e.g., [START_REF] Forrester | Blood parasites of raptors in Florida[END_REF][START_REF] Castellucci | Hematology and detection of hemoparasites in the American Kestrel (Falco sparverius) during summer nesting period[END_REF], the genera Haemoproteus and Leucocytozoon being the most frequently detected while Plasmodium commonly has lower prevalence [START_REF] Remple | Intracellular Hematozoa of raptors: a review and update[END_REF]). Overall, nestlings have been reported to exhibit relatively low prevalence of Haemosporida in numerous studies (e.g., [START_REF] Gutiérrez-López | Low prevalence of blood parasites in a long-distance migratory raptor: the importance of host habitat[END_REF][START_REF] Hanel | Blood parasites in Northern Goshawk (Accipiter gentilis) with an emphasis to Leucocytozoon toddi[END_REF]) even if notable exceptions exist [START_REF] Toyne | Blood parasites of nestling goshawks[END_REF]Ashford, 1997, Jeffries et al. 2015). Such differences of prevalences between juveniles and adults have been explained by the higher probability of adults contracting an infection in successive year [START_REF] Podmokła | Determinants of prevalence and intensity of infection with malaria parasites in the Blue Tit[END_REF]. For the Red Kite, available data does not allow to identify a difference of prevalence between adults and nestlings mainly because of the low number of positive individuals and unbalanced samples (Table 2).

The comparison of Haemosporida prevalence of Red Kite in Spain to those determined in other studies (Table 2) suggests that it could be higher in more northern areas such as France and potentially in Denmark, UK and Western Europe despite low sample sizes. The prevalence of haemosporidians in diurnal raptors has been explained by different factors. Some are related to populations such as the abundance of Diptera vectors within the environment [START_REF] Sol | Parasite mediated mortality and host immune response explain age-related differences in blood parasitism in birds[END_REF][START_REF] Sehgal | Manifold habitat effects on the prevalence and diversity of avian blood parasites[END_REF] while others rely more on species features, such as the world geographic range, breeding habitat and embryonic development [START_REF] Tella | Habitat, world geographic range, and embryonic development of hosts explain the prevalence of avian Hematozoa at small spatial and phylogenetic scales[END_REF]). In the case of the Red Kite, abundance of vectors has never been investigated directly in the previous studies published but it was assumed that the low prevalence measured in Spain could be due to the dryness of the environment, making it unsuitable for the main vectors of haemosporidian parasites [START_REF] Tella | Habitat, world geographic range, and embryonic development of hosts explain the prevalence of avian Hematozoa at small spatial and phylogenetic scales[END_REF]. In Franche-Comté, the climate is more wet (annual rainfall from 1000 to 1500 mm against 600 mm in Spain) and numerous ponds, lakes and streams composed the landscape of the studied areas, which could at least partly explain the relatively high prevalence of Leucocytozoon. Moreover, the world range of the Red Kite limited to Western Europe, the absence of inter-continental migration notably through Africa and the fact that it breeds in open-woodland habitat are all together factors that may explain relatively low Haemosporida prevalence in raptors and more generally animals [START_REF] Tella | Habitat, world geographic range, and embryonic development of hosts explain the prevalence of avian Hematozoa at small spatial and phylogenetic scales[END_REF][START_REF] Altizer | Animal migration and infectious disease risk[END_REF].

Lower body condition has been related to infection status (i.e., infected or not) by haemosporidians in both passerines [START_REF] Garvin | Blood parasites of nearctic-neotropical migrant passerine birds during spring trans-gulf migration: impact on host body condition[END_REF], Hatchwell et al. 2000) and raptors [START_REF] Dawson | Effects of hematozoan parasites on condition and return rates of American kestrels[END_REF]. We did not show any effect of infection of Red Kite by Leucocytozoon on body condition, which suggests that these blood parasites might have no or weak pathogenic consequences on nestlings and/or that the infection intensity is too low to impair condition. In the American kestrel Falco sparverius, [START_REF] Dawson | Effects of hematozoan parasites on condition and return rates of American kestrels[END_REF] found that a reduced condition of males was associated with higher intensity of infections by Haemoproteus using blood smear as a diagnostic tool. Moreover, [START_REF] Schoenle | An experimental test of the physiological consequences of avian malaria infection[END_REF] reported that measuring intensity infection rather than only presence/absence could be important when evaluating the role of infection in influencing hosts' behaviour, physiology or fitness. If the molecular methods we used are considered as relevant to diagnose infection status of Red Kite by blood parasites [START_REF] Martinsen | Morphological versus molecular identification of avian Haemosporida: An exploration of three species concepts[END_REF][START_REF] Krone | Haemosporidian blood parasites in European birds of prey and owls[END_REF], blood smears remain useful to measure the infection intensity in individuals and should be implemented complementarily to PCR screening in further studies.

Overall, both the low to medium prevalences reported in the different European populations studied and the absence of negative correlation with body condition we show suggest that hematozoans in Red Kite nestlings could be a limited threat for this life stage even if the detection of P. relictum in one chick questions about its occurrence at a larger scale and possible health effects on some individuals. In addition, recent findings showed that low-level chronic malaria infection reduced life span as well as the lifetime number and quality of offspring in a wild population of great reed warblers Acrocephalus arundinaceus, these delayed fitness effects of malaria being mediated by telomere degradation [START_REF] Asghar | Hidden costs of infection: Chronic malaria accelerates telomere degradation and senescence in wild birds[END_REF]. This stresses the need for further research on "subtle" hematozoan effect on individual in different avian groups including raptors. Moreover, the prevalence in adults has been only assessed in Spain, a country where it is particularly low, while breeder's survival is a critical parameter in the dynamic of raptor populations [START_REF] Saether | Avian life history variation and contribution of demographic traits to the population growth rate[END_REF].

Further investigation associating both qualitative and quantitative methods on a significant number of Red Kites and including adults, from other European areas, notably in Germany and France, would be needed to confirm whether blood pathogens can be considered as negligible threat for the Red Kite population.

Figure caption

0 [0] 0 [0] 2 Juvenile 2 0 [0] 0 [0] 0 [0] Adult 9 0 [0] 0 [0] 0 [0] nr 2 0 [0] 0 [0] 0 [0] 3 UK nr 1 0 [0] 0 [0] 0 [0] 4 Spain nr 46 0 [0] 0 [0] 0 [0] 5 Juvenile 164 4 [6] 0 [0] 0 [0] 6 Adult 104 2 [2] 0 [0] 0 [0] Denmark Nestling 17 18 [3] 0 [0] 0 [0] 7 Juvenile 1 0 [0] 0 [0] 0 [0] Adult 1 0 [0] 0 [0] 0 [0]

  et du Développement Durable (Plan National d'Actions pour la conservation du Milan Royal), by the Région Franche-Comté (BIOSCOPE project, convention 2016-0059) and by the Long Term Ecological Research (LTER) site Zone Atelier Arc Jurassien. Conflicts of interest/Competing interests: The authors declare that they have no interest or competing conflicts Ethical approval: All the manipulations were performed with authorizations no. 2015093-0007 and 2015093-0023 in the framework of the Red Kite national conservation plan. All the authors consent to participate and for publication Availability of data and material: The data presented in this article are referenced in the metadata platform dat@OSU proposed by the Observatoire des Sciences de l'Univers Terre Homme Environnement Temps Astronomie (OSU THETA) Code availability: NA Authors' contributions: Michael Coeurdassier designed the study. Christophe Morin and Michael Coeurdassier managed and carried out fieldwork and sampling. Eve Afonso and Nicolas Bassin conducted the pathogen analysis. Michael Coeurdassier, Eve Afonso, Thibaut Powolny and Nicolas Bassin analyzed and interpreted the data and wrote the manuscript. All co-authors commented on the manuscript and accepted the final version.

  linear model [body mass (in g) ~ wing length (in mm)], which corresponds to an assessment of the body mass following a statistical elimination (sensus Grafen and Hails 2002) of age differences. Positive residuals were related to chicks in relative good condition while negative residuals represented kites in poor condition. DNA extraction, PCR amplification, sequencing, and parasite detection DNA was extracted from blood samples using the DNeasy Blood & Tissue kit (Qiagen, France) and according to manufacturer's instructions (protocol "Purification of Total DNA from Animal Blood or Cells"). A nested PCR was conducted to amplify a fragment of cytochrome b gene (542 bp) from species of Plasmodium, Haemoproteus and Leucocytozoon as recommended by Pérez-Rodríguez et al. (2013). A first PCR was performed in total volumes of 12.5µL including 1 × Accustart II PCR ToughMix (Quanta Biosciences, USA), 0.6 mM of each primer (Plas1F/HaemNR3), and 22 ng of total genomic DNA. The PCR program consisted in an activation step of 3 min at 94°C, followed by 20 cycles of 30 sec at 94°C, 20 sec at 52°C, 30 sec at 72°C, and a final elongation of 5 min at 72°C. We then used PCR products from the first PCR to conduct the following reaction in a total volume of 25µL comprising 1 × Accustart II PCR ToughMix, 0.6 mM of each primer (3760F/HaemJR4) and 2.5µL of PCR products from the first PCR. The PCR program consisted in an activation step of 3 minutes at 94°C, followed by 35 cycles of 30 sec at 94°C, 20 sec at 52°C, 30 sec at 70°C, and a final elongation of 5 min at 70°C. PCR products were then separated and visualized using a Qiaxcel device (Qiagen, France) and a DNA High Resolution kit (Qiagen). A negative control (PCR grade water instead of DNA extract) was introduced every ten samples. PCR products which yielded an approximately 590 bp fragment (including primers) were confirmed by sequencing. The amplified products were purified using the QIAquick PCR Purification Kit (Qiagen) according to manufacturer's instructions. Direct sequencing of the purified PCR products was performed with an automated sequencer (Applied Biosystems 3130 Genetic Analyzer) using the primers employed for the second PCR reactions. For each PCR-positive sample, we aligned forward and reverse DNA sequences to obtain a consensus of 479 bp (Pérez-Rodríguez et al. 2013). A homology search of the sequences generated in this study was then performed by conducting an online search with the NCBI Basic Local Alignment Search Tool for nucleotides (BLASTN) in the GenBank database (http://blast.ncbi.nlm.nih.gov/Blast.cgi).Mixed infections were checked by examining electropherograms using the software BioEdit following the recommendations ofPérez-Tris and Bensch (2005). When double peaks were observed, we compared double peak patterns with the DNA sequences of parasite haplotypes already available on Genbank to identify which parasite haplotypes were involved in mixed infections[START_REF] Pérez-Rodríguez | Molecular characterization of haemosporidian parasites from kites of the genus Milvus (Aves: Accipitridae)[END_REF].Data analysisDifference of infection status (i.e., presence/absence of hematozoans) between sex was checked with a binomial generalized linear mixed effect model (GLMER) with 'Sex' as fixed factors and 'Nest' nested into 'Area' (Area/Nest) as random effects [corresponding full model: Status ~ Sex + (1|Area/Nest)]. The relationship between the presence of at least one positive chick in the brood ('Positive brood', binomial variable) and the brood size was tested with a binomial generalized linear mixed effect model (GLMER) with 'Brood size' as fixed factors and 'Area' as random effects [corresponding full model: Positive brood ~ Brood size + (1|Area)]. Relationship between body condition and infection status of the individual was checked with a linear mixed effect model (LMER) with 'Wing length', 'Sex' (in interaction or not with 'Wing length') and 'Status' as fixed factors and 'Nest' nested into 'Area' (Area/Nest) as random effects [corresponding full model: Body mass ~ Wing length + Sex + Status + (1|Area/Nest)]

Fig. 1

 1 Fig. 1 Body mass of the 47 Red Kite chicks studied standardised against wing

  (1981), 2: Krone et al. (2001), 3: Krone et al. (2008), 4: Peirce et al. (1983), 5: Tella et al. (1999), 6: Pérez-Rodriguez et al. (2013), 7: Laursen (2008).

Table 1

 1 Linear mixed-effect models tested to explain body mass as response. The best model selected based on corrected Akaike information criteria is reported in bold (LL: maximized log-likelihood; K: number of estimated parameters; AICc: Corrected Akaike's information criterion; w_ic: Akaike weight)

	Model tested (Body mass ~ )	LL	K	AICc	ΔAICc	w_ic
	Wing length + Sex + (1|Area/Nest)	-258.9	6	532.0	0.00	0.39
	Wing length + Sex * Status + (1|Area/Nest)	-256.5 8	532.9	0.89	0.25
	Wing length + Sex + Status + (1|Area/Nest) -258.1 7	533.0	1.03	0.23
	Wing length * Sex + (1|Area/Nest)	-258.7	7	534.2	2.20	0.13
	Wing length + (1|Area/Nest)	-266.1 5	543.6	11.57	0.00

Table 2

 2 Summary of the data available on Red Kite infection by the haemosporidian parasites Leucocytozoon spp. (Leuco.), Haemoproteus spp. (Haemop.) and Plasmodium spp. (Plasmo.). (n = total number of screened individuals; nr = not reported; ns = not screened)

	Area/country Age	n	Prevalence in % [n of infected kites]	Reference
				Leuco.	Haemop.	Plasmod.
	Western	nr	4	25 [1]	25 [1]	0 [0]	1
	Europe
	Germany	Nestling	13	15 [2]

We warmly thank the three anonymous referees for their contribution to improve this article. This study was financially supported by the French Ministère de l'Ecologie et du Développement Durable (Plan National d'Actions pour la conservation du Milan Royal), by the Région Franche-Comté (BIOSCOPE project, convention 2016-0059) and by the LTER site Zone Atelier Arc Jurassien. The molecular analyses were carried out on the PEA²t platform. The experiments comply with the current laws of France when they were performed.


Valentin Goussard

Francois Duprat

Vincent Gerbaud

Jean-Luc Ploix

Geŕard Dreyfus

Veŕonique Nardello-Rataj

Jean-Marie Aubry

Predicting the Surface Tension of Liquids: Comparison of Four Modeling Approaches and Application to Cosmetic Oils

The efficiency of four modeling approaches, namely, group contributions, corresponding-states principle, σ-moment-based neural networks, and graph machines, are compared for the estimation of the surface tension (ST) of 269 pure liquid compounds at 25 °C from their molecular structure. This study focuses on liquids containing only carbon, oxygen, hydrogen, or silicon atoms since our purpose is to predict the surface tension of cosmetic oils. Neural network estimations are performed from σ-moment descriptors as defined in the COSMO-RS model, while methods based on group contributions, corresponding-states principle, and graph machines use 2D molecular information (SMILES codes). The graph machine approach provides the best results, estimating the surface tensions of 23 cosmetic oils, such as hemisqualane, isopropyl myristate, or decamethylcyclopentasiloxane (D5), with accuracy better than 1 mN•m -1 . A demonstration of the graph machine model using the recent Docker technology is available for download in the Supporting Information.

INTRODUCTION

Surface tension (ST) is the force existing at a liquid/gas interface, which tends to minimize the surface area; it is caused by asymmetries in the intermolecular forces between molecules located at the interface. It plays an important role in a number of processes where a liquid/gas interface is present; for instance, it drives the shape of small liquid drops and bubbles or the wetting of a solid surface by a liquid. [START_REF] Delgado | A Molecular Structure Based Model for Predicting Surface Tension of Organic Compounds[END_REF] In cosmetics, the spreadability of oils on the human skin reflects both the sensory qualities and the efficacy of the product; 2,3 Parente et al. reported that the spreadability of oils and the film forming properties could be partly correlated to the surface tension of the cosmetic oil: 2 the lower the surface tension, the more spreadable the oil.

Surface tensions can be measured by numerous methods (e.g., Wilhelmy plate, Du Nouÿ ring, maximum bubble pressure, dropweight, or hanging-drop) accounting for the large amount of experimental data available in the literature. [START_REF] Gharagheizi | Handling a Very Large Data Set for Determination of Surface Tension of Chemical Compounds Using Quantitative Structure-Property Relationship Strategy[END_REF] However, for a given compound, significantly different surface tension values can be found, depending on the method used and/or the purity of the liquid. Furthermore, when virtual liquids are generated in silico, [START_REF] Gharagheizi | Use of Artificial Neural Network-Group Contribution Method to Determine Surface Tension of Pure Compounds[END_REF] the experimental determination of surface tension is not possible. In this case, it is crucial to predict reliably the surface tension values of candidate compounds in order to select the most promising ones, which should be synthesized first. Unfortunately, predicting this property is not straightforward because of the complexity of the phenomenon, i.e., the necessity to break bonds in order to bring a molecule from the bulk phase to the surface. 6 Various prediction methods for surface tension values were described previously by several authors. [START_REF] Delgado | A Molecular Structure Based Model for Predicting Surface Tension of Organic Compounds[END_REF][6][START_REF] Conte | Combined Group-Contribution and Atom Connectivity Index-Based Methods for Estimation of Surface Tension and Viscosity[END_REF][8][START_REF] Kauffman | Prediction of Surface Tension, Viscosity, and Thermal Conductivity for Common Organic Solvents Using Quantitative Structure-Property Relationships[END_REF][START_REF] Freitas | A Linear Free Energy Analysis of the Surface Tension of Organic Liquids[END_REF][START_REF] Albahri | Modeling of Pure Compounds Surface Tension Using QSPR[END_REF] In this paper, four modeling approaches, namely, group contributions, corresponding-states principle, σ-moment-based neural networks, and graph machines, are applied to the prediction of the surface tension of liquid compounds. These methods are described below.

Molecular dynamic simulations were first used in 1976, when Rao et al. highlighted the possibility of computing the surface tension of a liquid film from statistical thermodynamic relations. [START_REF] Rao | Surface Structure of a Liquid Film[END_REF] In 2016, Ghoufi et al. calculated the surface tension by modeling the liquid/gas interface. [START_REF] Ghoufi | Computer Modelling of the Surface Tension of the Gas-liquid and Liquid-liquid Interface[END_REF] Density functional theory has also been used successfully for the surface tension prediction of 18 pure nonpolar fluids. [START_REF] Fu | Prediction of Surface Tension for Pure Non-Polar Fluids Based on Density Functional Theory[END_REF] Nevertheless, models based on molecular dynamics are computer demanding and not suitable for structures with more than a few tens of atoms such as some we are interested in.

Introduced by Brock and Bird in 1955, [START_REF] Brock | Surface Tension and the Principle of Corresponding States[END_REF] the correspondingstates principle allows surface tension estimations from the critical constants P c , V c , and T c . Curl and Pitzer improved this model, which is restricted to short and weakly polar molecules by introducing an acentric shape parameter. [START_REF] Curl | Volumetric and Thermodynamic Properties of Fluids-enthalpy, Free Energy, and Entropy[END_REF] Escobedo et al. used this model for the surface tension prediction of 94 compounds, [START_REF] Escobedo | Surface Tension Prediction for Pure Fluids[END_REF] and Li et al. presented a new method based on a combination of corresponding-states and group contribution models. [START_REF] Li | Estimations of Surface Tensions at Different Temperatures by a Corresponding-States Group-Contribution Method[END_REF] Group contribution-based (GC) methods are popular approaches to predict a wide range of physicochemical properties. [START_REF] Marrero | Group-Contribution Based Estimation of Pure Component Properties[END_REF][START_REF] Constantinou | New Group Contribution Method for Estimating Properties of Pure Compounds[END_REF] In such methods, the molecular structure of any organic compound is viewed as an assembly of functional groups with specific property attributes. The properties of the studied molecule are computed as a function of the contributions of the different groups, with respect to their occurrence (firstorder contribution), with respect to adjacent groups (second order), and with respect to other molecular attributes (third order). Regarding the prediction of surface tension, a model was introduced by Egemen et al. in 2000 to predict the surface tension of 44 organic compounds, [START_REF] Egemen | Predicting Surface Tension of Liquid Organic Solvents[END_REF] and an extended model was proposed by Conte et al. [START_REF] Conte | Combined Group-Contribution and Atom Connectivity Index-Based Methods for Estimation of Surface Tension and Viscosity[END_REF] Models based on quantum chemistry are also effective to predict a wide range of physicochemical properties of organic compounds. In particular, the COSMO-RS method (Conductorlike Screening MOdel for Real Solvents), introduced by Klamt in 1995, is based on unimolecular quantum chemical calculations combined with exact statistical thermodynamics. [START_REF] Klamt | From Quantum Chemistry to Fluid Phase Thermodynamics and Drug Design[END_REF][START_REF] Klamt | Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena[END_REF][START_REF] Klamt | Refinement and Parametrization of COSMO-RS[END_REF] It provides the necessary information for predicting chemical potentials in liquids, from which one can compute many other thermodynamic equilibrium properties such as activity coefficients, solubility, partition coefficients, vapor pressure, and free energy of solvation, [START_REF] Klamt | A Novel View to Physiological Solvation and Partition Questions[END_REF] without requiring any experimental input. For instance, the interfacial tension between two liquids, i.e., the free energy per unit surface area that is required for maintaining an interface between two condensed phases, can be calculated using COSMO-RS. [START_REF] Andersson | First-Principles Prediction of Liquid/Liquid Interfacial Tension[END_REF] However, the prediction of surface tensions remains a challenge, yet unsolved by COSMO-RS (Klamt personal communication). Interestingly, Kondor et al. proposed to predict surface tensions at different temperatures by a regression method, using the so-called Klamt set of COSMO-RS σ-moments. This set of five descriptors is well known to encode as much chemical information as the complete set of 15 COSMO-RS σ-moments. The five σ-moments calculated from the σ-profiles of the molecules under study [START_REF] Klamt | From Quantum Chemistry to Fluid Phase Thermodynamics and Drug Design[END_REF][START_REF] Klamt | A Novel View to Physiological Solvation and Partition Questions[END_REF][START_REF] Lukowicz | Rationalization and Prediction of the Equivalent Alkane Carbon Number (EACN) of Polar Hydrocarbon Oils with COSMO-RS σ-Moments[END_REF] and the temperature were selected as inputs of a nonlinear Quantitative Structure-Property Relationship (QSPR) model, 6 based on neural networks.

QSPR models are regression or classification models that predict a property of a molecule from selected descriptors thereof; they were used successfully for the prediction of surface tensions of a wide range of organic compounds. [START_REF] Delgado | A Molecular Structure Based Model for Predicting Surface Tension of Organic Compounds[END_REF] [START_REF] Kauffman | Prediction of Surface Tension, Viscosity, and Thermal Conductivity for Common Organic Solvents Using Quantitative Structure-Property Relationships[END_REF] Freitas et al. decreased the number of descriptors of the model to six to successfully predict the surface tensions of 299 compounds. [START_REF] Freitas | A Linear Free Energy Analysis of the Surface Tension of Organic Liquids[END_REF] Delgado et al. designed a model with six descriptors corresponding to the different intermolecular interactions in the bulk phase and calculated from the molecular structures 1 for the prediction of the surface tensions of 320 compounds.

Goulon et al. developed graph machines as an alternative QSAR/QSPR method. [START_REF] Goulon | Predicting Activities without Computing Descriptors: Graph Machines for QSAR[END_REF] Graph machines are regression or classification models that estimate properties of molecules directly from the topological information supplied by their SMILES (Simplified Molecular Input Line Entry Specification) codes. Therefore, they do not require any descriptor. In these models, molecules are described as graphs derived from their 2Dstructure, and the parametrized functions that compute the estimation of the property or activity of interest reflect the compound molecular structures. As usual in regression or classification models, the parameters of the graph machines are computed during training from examples present in a database of experimental values of the property or activity of interest. The main benefit of the method is the absence of descriptors, with the SMILES codes being the only required information.

In the present work, surface tension estimations are performed by four different methods, falling into two families: (i) group contribution models and corresponding-states principle models (Section 2.1) and (ii) neural network regression and graph machine regression (Section 2.2). Neural network regression is performed from COSMO-RS σ-moments, whose computation is described in Section 2.3. For the first family, methods are taken from the literature. For the two methods of the second family, predictive models are designed and trained from a database of surface tensions, at 25 °C, of 269 molecules belonging to several chemical families and containing only carbon, hydrogen, oxygen, and silicon atoms, namely, alkanes, ethers, esters, ketones, carbonates, acids, alcohols, silanes, and siloxanes, with straight, branched, or cyclic chains. For all four methods, the predictive ability of the resulting models are assessed by the estimation of the surface tensions of 23 cosmetic oils that are not present in the training data set.

EXPERIMENTAL AND THEORETICAL METHODS

2.1. Group Contribution and Corresponding-States Based Models. Surface tension estimations using methods based both on group contribution (GC) or on the corresponding-states principle (CSP) were carried out using the InBioSynSolv (IBSS) property library developed by Heintz et al. [START_REF] Heintz | Computer Aided Product Design Tool for Sustainable Product Development[END_REF] IBSS software is intended to design new molecules satisfying targeted properties and to predict many physicochemical properties by means of methods encoded in a property model library. [START_REF] Conte | Combined Group-Contribution and Atom Connectivity Index-Based Methods for Estimation of Surface Tension and Viscosity[END_REF][START_REF] Heintz | Computer Aided Product Design Tool for Sustainable Product Development[END_REF] Among the methods available in IBSS for the prediction of surface tensions, the Conte et al. GC method [START_REF] Conte | Combined Group-Contribution and Atom Connectivity Index-Based Methods for Estimation of Surface Tension and Viscosity[END_REF] was developed by splitting 402 molecules into chemical groups and estimating the parameters of the following relation:

∑ ∑ ∑ = + + ST NC MD OE i i i j j j k k k ( 1
)
where ST is the surface tension, C i , D j , and E k are the first-, second-, and third-order contributions estimated by successive linear regressions, and N i , M j , and O k are the occurrences of the molecular groups. The reported standard deviation of the estimation error (i.e., of the difference between the estimated and measured values of the surface tension) is 1.47 mN•m -1 , and the mean absolute estimation error is 1.05 mN•m -1 over the 402 data points of the database used for estimating the parameters C i , D j , and E k , in a range between 15 and 52 mN•m -1 . [START_REF] Conte | Combined Group-Contribution and Atom Connectivity Index-Based Methods for Estimation of Surface Tension and Viscosity[END_REF][START_REF] Hukkerikar | Group-Contribution + (GC + ) Based Estimation of Properties of Pure Components: Improved Property Estimation and Uncertainty Analysis[END_REF] Once the contribution values are estimated, the GC method is intended to predict the property value of any molecular structure. Difficulties can appear when using such methods, e.g., a nonunique decomposition of some structures or compounds with fragments for which no group contributions have been defined yet. This is the case for siloxane derivatives, which are at the core of the present study.

Other estimation methods have been developed, based on the principle of corresponding states. [START_REF] Poling | The Properties of Gases and Liquids[END_REF] This principle states that some dimensionless properties will follow the same relation over dimensionless variables of state and other dimensionless quantities for all compounds, typically, the reduced temperature (T r = T/T c ) and reduced pressure (P r = P/P c ), where T c and P c are the critical temperature and pressure. The related CSP models are the so-called two-parameter corresponding-states principle relation. For surface tension, an additional parameter, Pitzer's dimensionless acentric factor ω, was introduced to better match the CSP relation and the experimental data, leading to the following relation: [START_REF] Curl | Volumetric and Thermodynamic Properties of Fluids-enthalpy, Free Energy, and Entropy[END_REF][START_REF] Pitzer | [END_REF] (3)

ω ω ω = + + - - ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ST P T T 1.
ω ω ω ω = + - - - ST ST ST ST ( ) r r r r (1) (1)
(

(

2) (1) (4)
where

= - ST T 40.520(1 ) r r (1) 1.287 (5)
In relations 2 and 3, with the critical pressure and temperature expressed in bar and Kelvin, respectively, ST is obtained in mN• m -1 . The authors claimed a five percent accuracy for nonpolar fluids and reported larger deviations for more polar substances.

2.2. Neural Network and Graph Machine Modeling. 2.2.1. Database Construction. The data set required for the training of neural networks and graph machines was carefully built with the following requirements in mind: (i) to obtain a representative sample of the surface tension variation with the chain length for compounds belonging to homologous families, (ii) to select molecules that contain the constitutional atoms usually found in the cosmetic oils, i.e., carbon, hydrogen, oxygen, and silicon, (iii) to select, for each molecule, at least two consistent measured values for the surface tension provided by different literature references, and (iv) to include only molecules that are in the liquid state at 25 °C under atmospheric pressure. A data set of 269 molecules (hereinafter termed "complete set") with a measured surface tension ranging from 12 to 48 mN•m -1 was built from data extracted from reliable sources [START_REF] Jasper | The Surface Tension of Pure Liquid Compounds[END_REF][START_REF] Le Neindre | Tensions superficielles des composeś organiques[END_REF][START_REF] Yaws | Chemical Properties Handbook: Physical, Thermodynamic, Environmental, Transport, Safety, and Health Related Properties for Organic and Inorganic Chemicals[END_REF][START_REF] Wohlfarth | Numerical Data and Functional Relationships in Science and Technology[END_REF][START_REF] Haynes | CRC Handbook of Chemistry and Physics[END_REF] and retrieved from technical datasheets. [START_REF]Comparative Solvents Data[END_REF][START_REF][END_REF] When measured data were not available at 25 °C, but for several other temperatures nearby, a simple linear equation was used to interpolate the surface tension at 25 °C.

Model Selection.

Model selection is a crucial step in the design of machine-learning based models. Its purpose is to find, given the available data, the appropriate model complexity that provides the best generalization, i.e., that provides the most accurate estimations on data that are not present in the training set. The available data set of 269 molecules was partitioned into a training/validation set of 244 molecules and a test set of 25 molecules for performance assessment. The molecules of the test set were chosen so that (i) the distributions of molecules among the chemical families considered were very similar in the complete and test sets and that (ii) the distribution of the surface tension values was as uniform as possible on the range of measured values. The distribution of the molecules among chemical families is displayed in Figure 1.

Model selection was performed by computing the virtual leave-one-out score (defined in Section 3.2) of each model, at the end of the training process performed with the training/ validation set. Further details are provided in Sections 2.2.3 and 2.2.4.

After model selection, the performance of the selected neural network model and the performance of the selected graph machine model were assessed on the 25-molecule test set. As a final step, the surface tensions of 23 cosmetic oils were measured as described in Section 2.4 and predicted by the neural network and the graph machine models as described in Section 3.3.

Neural Networks.

A neural network is a nonlinear parametrized function. [START_REF] Dreyfus | Neural Networks: Methodology and Applications[END_REF] In the framework of nonlinear regression for estimating a quantity of interest from a training set of N T elements, the neural network parameters are estimated from available data during a training phase, which consists of minimizing the sum of squared training errors J(θ) with respect to the parameters

∑ θ θ = - = J y g ( ) ( ( )) i N i i 1 2 T (7)
where y i is the measured value of the quantity of interest for the ith element of the training set, and g i (θ) is the value of the quantity of interest estimated by the neural network. In the present case, the quantity of interest is the surface tension of liquids, hence y i is the measured value of the surface tension of molecule i of the database, g i (θ) is the estimated value of the surface tension of molecule i, and N T = 244 as explained in the previous section. In the present study, the minimization of J(θ) was performed by the Levenberg-Marquardt algorithm.

For nonlinear regression by neural networks, the usual form of g i (θ) is a linear combination of nonlinear functions called "hidden neurons". Each hidden neuron is a continuous, nonlinear function (typically a tanh function) of a linear combination of the variables. Such a neural network is called a multilayer Perceptron (MLP). It is particularly suitable for nonlinear regression because it is a universal approximator: any continuous nonlinear function can be approximated with arbitrary accuracy by an MLP having a finite number of hidden neurons. An MLP with zero hidden neurons is a linear model, hence is trained by traditional least-squares fitting.

As usual in nonlinear regression, the generalization ability of an MLP model depends on its complexity, which depends basically on the number of hidden neurons. A neural network with too many neurons approximates the training data very accurately but is unable to generalize, i.e., to predict accurately data that are not present in the training set; a neural network with too small a number of neurons is unable to learn the training data, hence to generalize (this phenomenon is known as the "bias-variance dilemma"). Therefore, the design of an MLP requires a model selection phase. [START_REF]Comparative Solvents Data[END_REF] In the present study, MLP selection was performed by virtual leave-one-out, as described in Section 3.2.

2.2.4. Graph Machines. As mentioned previously, a graph machine is intended to perform classification or regression on graphs; in QSAR/QSPR applications, the graphs under consideration are the 2D structures of the molecules. In the present case, a graph machine provides an estimate of the surface tension; the latter quantity being continuous, the task considered here is a regression task.

The design of a graph-machine based model includes the following steps:

• Construction of the 2D-graph of the molecule from its SMILES representation: Each node of the graph is a non-H atom, and each edge is a chemical bond of the graph.

Each node has at least two labels: its nature and its degree (the number of chemical bonds that bind it to its adjacent non-H atoms). For molecules that contain cycles, hence are represented by a cyclic graph, edges are deleted in order to form an acyclic graph in which every path of the graph has its end at a specific node called "output node". [START_REF] Dioury | QSPR Prediction of the Stability Constants of Gadolinium(III) Complexes for Magnetic Resonance Imaging[END_REF] and references therein. The bias-variance dilemma applies to graph machines as well as to any machine learning method, so an appropriate model selection procedure must be performed. The selection method for graph-machine based models is described in Section 3.2.

2.3. COSMO-RS Computing Procedure. COSMO-RS uses a combination of quantum chemistry and statistical thermody-namics in order to compute the chemical potential μ of a solute in a liquid phase, which can then be converted into physicochemical properties. [START_REF] Klamt | From Quantum Chemistry to Fluid Phase Thermodynamics and Drug Design[END_REF] As mentioned in the Introduction, no model for the estimation of surface tension is available from COSMO-RS yet. An alternative approach consists in using descriptors provided by COSMOtherm: the so-called COSMO-RS σ-moments.

Because of polar covalent bonds, molecules bear a charge density σ on the so-called σ-surface, which is a slightly inflated van der Waals surface of the compound. The full 3D information about the charge density repartition on the σ-surface can be reduced to a curve p X (σ), the σ-profile of the molecule X, which is a smoothed histogram expressing how much of the σ-surface lies in the polarity interval [σdσ/2, σ + dσ/2]. [START_REF] Klamt | A Novel View to Physiological Solvation and Partition Questions[END_REF] Klamt suggested an expression of the partition coefficient K of a solute X between two phases as a development on the σmoments (eq 8). [START_REF] Klamt | From Quantum Chemistry to Fluid Phase Thermodynamics and Drug Design[END_REF] Usually, a development to m = 6 is sufficient to estimate satisfactorily the partition coefficient K.

∑ = + + = RT K c M c M cM ln acc acc X don don X i m i i X 0 (8)
The σ-moments M i X are calculated from the σ-profiles p X (σ) of the studied compound X, according to eqs 9-11:

∫ σ σ σ σ = - σ + +∞ M p d ( )( ) acc X X HB HB (9) ∫ σ σ σ σ = -- σ -∞ - M p d ( )( ) don X X HB HB (10) ∫ σ σ σ = -∞ -∞ M p d ( ) i X X i ( 11
)
Among the 15 σ-moments, the first ones have a simple physical meaning. The zero-order σ-moment M 0 X is the whole surface area of the molecule, expressed in Å 2 . The first-order one M 1 X is the polarization charge of this surface, expressed in e. For uncharged molecules, this moment is equal to zero. The second-order σ- , expressed in e 2 •Å -2 , is the polarity of the studied molecule. [START_REF] Lukowicz | Rationalization and Prediction of the Equivalent Alkane Carbon Number (EACN) of Polar Hydrocarbon Oils with COSMO-RS σ-Moments[END_REF] The third-order M 3 X represents the asymmetry of the σ-profile p X (σ). The other σ-moments M 4 X , M 5 X , and M 6 X have no known physical meaning. Finally, M acc X and M don X , expressed in e, are the "hydrogen-bonding" σ-moments. They represent the ability of the molecule to interact with hydrogen-bond acceptors and donors, respectively. They have a nonzero value when the σprofile of a molecule outstrips the range [-σ HB , +σ HB ], where σ HB is the hydrogen-bond cutoff equal to 0.01 e•Å -2 .

The 14not 15 since M 1 X is always equal to zero for the uncharged molecules of our studycalculated σ-moments values for the complete and test sets of molecules are given in Tables S3 andS4 of the Supporting Information.

2.4. Surface Tension Measurements for Cosmetic Oils. To assess the accuracy of the estimations of the developed model, a test data set of 23 cosmetic oils that are not present in the complete set was built. The surface tensions of these cosmetic oils were determined experimentally, using two different methods: measurements were both performed (i) with the Kruss K100 tensiometer which measures the force occurring when a probe (ring, plate or rod) is completely wetted with a liquid and (ii) by the pendant drop method associated with the Kruss DSA10 tensiometer. It has thus been verified on a few molecules that values obtained for surface tensions are consistent from one method to the other. The values measured with the K100 tensiometer at 23.5 °C are presented in Table 3. The rod is immersed in each liquid studied, with an immersion rate of 10 mm•min -1 and an immersion depth of 2 mm. The force acting on the rod is recorded and allows the calculation of the surface tension of the liquid. Before each series of measurements, calibration with pure water is carried out. Additionally, surface tensions of four liquids belonging to the training/validation set, whose values (reported in the literature) can be found in Table S1 of the Supporting Information, were measured for the following results (all values with a standard deviation of 0.2 mN• m -1 ): n-undecane ST = 23.1 mN•m -1 , n-pentadecane ST = 26.9 mN•m -1 , triacetine ST = 35.1 mN•m -1 , and decamethylcyclopentasiloxane ST = 18.0 mN•m -1 . These values are consistent with those of the literature. During a liquid measurement, the surface tension values were sampled every 3 s until the standard deviation was below 0.1 mN•m -1 . Three measurements were taken in a row for each sample. The means and standard deviations of the measurements performed on the 23 compounds are reported in Table S2 of the Supporting Information. No temperature correction was made to the previous measured values since the observed deviations between 0.1 and 0.7 mN•m -1 were larger than the temperature corrections, expected to be below 0.2 mN•m -1 .

RESULTS AND DISCUSSION

3.1. Dependence of Surface Tension on Chemical Function and Chain Length. The variation of the surface tensions of several chemical families (alkanes, ethers, esters, ketones, carbonates, acids, alcohols, silanes, and siloxanes) with the number of carbon atoms is plotted in Figure 3, gathering experimental values from at least two independent sources. For silicon-containing compounds, the variation of the surface tension is plotted against the sum of the number of carbon and silicon atoms. For instance, decamethylcyclopentasiloxane (C 10 H 30 O 5 Si 5 ) is plotted at abscissa 15 on the X-axis of Figure 3.

For all chemical families, the surface tension increases with the number of carbon or (carbon+silicon) atoms. Surface tension is closely related to the intermolecular forces: the stronger these forces, the more tightly bound the molecules, hence the higher the surface tension. Because long chains are more likely to form more intermolecular bonds than small ones do, the surface tension increases with chain length. Silicones exhibit drastically low surface tension. This is due to their inorganic polar backbone on which apolar organic groups are grafted, able to rotate freely around the Si-O bonds. These apolar groups shield the Si-O backbone, preventing the dipolar interactions thereby decreasing the surface tension. Apart from silicones, alkanes and ethers have the lowest surface tensions because only London interactions take place between CH 2 /CH 3 groups. Compounds containing polar ester groups have higher surface tensions because of the Keesom and Debye interactions. Finally, alcohols have higher surface tension values because these molecules carry a free hydroxyl group that is involved in strong intermolecular hydrogen bonds. As explained in the Introduction, a low surface tension results in high spreadability. That is why silicones are used in cosmetic formulations. Figure 3 gives an indication on how chemical functions affect the surface tension of linear molecules, showing that alkanes and ethers are the carbon compounds whose surface tensions are closest to those of silicones. It also focuses on the effect of chain length on the surface tension for series of linear molecules (except cyclic silicones), highlighting the increase in intermolecular forces with the chain length, responsible for an increase in surface tension. Obviously, this property does not depend only on those two factors, but also on topological modifications such as cyclization or branching of compounds. Thus, results presented in the following sections of this study also deal with branched and cyclic compounds.

3.2. Modeling Surface Tensions of the Liquids of the Complete Set. 3.2.1. Surface Tension Estimation from σ-Moments by Neural Networks (Multilayer Perceptrons). In the present section, we discuss the design and test of a multilayer Perceptron for the estimation of the surface tension of the molecules of the complete data set described in Section 2.2.1, comprising 269 compounds that carry diverse functional groups (ethers, ketones, acids, carbonates, esters, silanes, siloxanes, alkanes, and alcohols) and that possess straight or branched chains, cycles, etc. The compounds contain only carbon, hydrogen, oxygen, and silicon atoms because those are the only atoms present in the backbone of cosmetic oils. The database was partitioned into a training/validation set of 244 compounds and a test set of 25 compounds as described in Section 2.2.3. In the spirit of Kondor et al., 6 the estimation was performed by a multilayer Perceptron from five COSMO-RS σmoments: M 0 , M 2 , M 3 , M acc , and M don . [Neural network computations were performed with a homemade software package based on the NeuroOne kernel, a former product of Netral S.A.] Feature selection by the random probe method [START_REF] Dreyfus | Neural Networks: Methodology and Applications[END_REF] confirmed the fact that for our data M 0 , M 2 , and M 3 are the most relevant for ST estimation.

Model selection was performed by training models of increasing complexity (number of hidden neurons). [START_REF]Comparative Solvents Data[END_REF] For each model complexity, 100 models were trained with different random initial parameter values. The Root Mean Square Training Error (RMSTE) and the Virtual Leave-One-Out (VLOO) score were computed:

θ = RMSTE J N ( ) m T ( 12
)
where J(θ) is the sum of square errors (eq 7), θ m is the vector of parameters of the model after completion of training, and N T = 244 is the number of molecules in the training/validation set:

∑ θ = - = - VLOO score N y g 1 [ ( )] T i N i i 1 m i 2 T (13)
where y i is the surface tension value measured for molecule i, and

g i (θ m -i
) is the surface tension value for molecule i that would have been estimated by the model if molecule i had been left out the training set (i.e., the model had been trained on all molecules of the training/validation set except molecule i). The virtual leaveone-out score provides an unbiased estimation of the generalization ability of the model. [START_REF] Monari | Local Overfitting Control via Leverages[END_REF] The 10 models (out of 100) having the smallest VLOO scores were stored.

The third row of Table 1 shows the values of the RMSTE on the 244 molecules of the training/validation set for the best trained models (models with minimum RMSTE) of increasing complexity. As expected, RMSTE decreases monotonously with increasing complexity; by contrast, the mean value of the VLOO scores of the 10 stored models (fourth row) first decreases, reaches a minimum for the model with 78 parameters (multilayer Perceptron with 11 hidden neurons), and subsequently increases. This illustrates the bias-variance dilemma explained in Section 2.2.3. The VLOO score difference between the 71-and 78-parameter models is very small, and the recommended practice in such a situation is to select the most parsimonious (least complex) model. Therefore, the multilayer Perceptron with 10 hidden neurons was selected for testing. The results of the models with less than four hidden neurons are not displayed, as they have much higher values.

For the selected modelwith 10 hidden neuronsthe VLOO ST estimates for molecules of the training/validation set and the ST estimates for the molecules of the test set are plotted against their measured values in Figure 4. The RMSE computed for the training/validation and test sets are equal to 1.91 and 3.60 mN•m -1 , respectively, and the determination coefficients R 2 , displayed for both sets, are equal to 0.89 and 0.76, respectively. The surface tension is correctly estimated around the center of the range but is clearly inadequate for the outer examples, which impairs the R 2 value for the test set where outliers, being at the ends of the range, have large leverages.

As mentioned in the Introduction, Kondor et al. were the first to use the COSMO-RS σ-moments for nonlinear regression by neural networks. To build their model, Kondor et al. used a database of 1500 surface tensions recorded for 188 molecules at different temperatures. 6 The initial set was partitioned into a training/validation set of 1275 examples and a test set of 225 examples. The results, obtained with a neural network whose variables were the COSMO-RS σ-moments mentioned above, were fair for both sets. However, the method used for partitioning the set remains unclear since there are more examples in the test set than the number of initial molecules. Therefore, some compounds are inevitably in both sets at least once, though at different temperatures. This is clearly identified for CO 2 , with three data points belonging to the training set and one to the test set. A solid assessment of the prediction ability of Kondor's model would request a test with fresh examples not pertaining to any of the previous sets and possessing known surface tension values at various temperatures. To the best of our knowledge, this is yet to be done.

Surface Tension Estimation from SMILES by Graph

Machines. In the present section, we discuss the design and test of graph machines (described in Section 2.2.4) for the estimation of the surface tension of the molecules of the complete data set described in Section 2.2.1, comprising 269 compounds, and partitioned into a training/validation set of 244 compounds and a test set of 25 compounds.

The model design methodology for graph machines is the same as described in Section 3.2.1 for multilayer Perceptrons: graph machines of increasing complexity (number of hidden neurons of the node function) were trained. For each model complexity, 100 models were trained with different random initial parameter values. RMSTE (eq 12) and the VLOO scores (eq 13) were computed, and the 10 models having the smallest VLOO scores were stored.

The third row of Table 2 shows the RMSTE values of the besttrained models (models with smallest RMSTE) with different complexities for the 244 molecules of the training/validation set. As expected, the RMSTE decreases with increasing complexity, while the mean value of the VLOO scores of the 10 stored models (fourth row) goes through a minimum (model with 95 parameters, i.e., node function with seven hidden neurons) and subsequently increases. Therefore, graph machines with seven hidden neurons were selected for further testing.

As an illustration, the VLOO ST estimates for molecules of the training/validation set and the ST estimates for the molecules of the test set are plotted against their measured values in Figure 5. The RMSE computed for the training/validation and test sets are both equal to 0.77 mN•m -1 , and the determination coefficients R 2 , displayed for both sets, are above 0.98. This demonstrates that the VLOO score on the training/validation set provides an accurate assessment of the generalization ability of the model. It also shows that the accuracy of the graph machine estimations is much better than the accuracy of MLP estimations from the σmoments. Indeed, the estimation error of graph machines is on the order of the experimental uncertainty on the available data.

Tables gathering the results with both group contribution and corresponding-states based methods, and respective RMSE, are given in the Supporting Information.

3.3. Comparison of the Four Approaches on the 23 Cosmetic Oils Test Set. For further testing, the cosmetic oils set (23 molecules) was investigated. Most molecules display polar functions with oxygen. As such we may expect large deviation in the prediction with the Zuo and Stenby corresponding-states method that was developed over a database containing only nonpolar molecules. To take advantage of the maximum information on all the available data, the complete set of 269 molecules was used as a training/validation set. For graph machines, 250 models with node functions having seven hidden neurons were trained from that set. The 25 models that had the smallest VLOO scores were applied to the cosmetic oils set, and the mean of their estimates of the surface tension was taken as the final estimate. This resulted in a RMS estimation error of 0.9 mN• m -1 . For the MLP estimation from σ-moments, 250 models having 10 hidden neurons were trained from the same set. The 25 models that had the smallest VLOO scores were also applied to the cosmetic oils set, and the mean of their estimates of the surface tension was taken as the final estimate. The results are compared with those obtained with the other two methods, group contribution, and corresponding-states method. Deviations from experimentally measured values of surface tension are gathered in Table 3. Because of the absence of contribution for siloxane groups in Conte's group contribution method, the surface tension for the dodecamethylcyclohexasiloxane (D6) cannot be estimated.

A large deviation is observed between the experimental data and the values predicted by the GC method. We notice that the molecules considered in Table 3 have a larger molecular weight than compounds included in the database used for the predictions. Conte's original paper did not give a list of data set molecules but provided several examples. [START_REF] Conte | Combined Group-Contribution and Atom Connectivity Index-Based Methods for Estimation of Surface Tension and Viscosity[END_REF] Not a single one concerned such large molecules. This hints at the fact that the database set used to develop the GC model does not sample the space of molecules to which the test set we are interested in belongs. Thus, in the following, surface tensions estimated with the corresponding-states methods described by Zuo and Stenby are compared to results obtained with the neural networks and graph machines. Figure 6 presents the comparison between the experimentally measured surface tensions of the 23 cosmetic oils of the test set and the predicted surface tensions with corresponding-states based method, neural networks using the COSMO-RS σ-moments as descriptors, and graph machines.

The method of corresponding states advocated by Zuo and Stenby performs poorly as was expected. For all but for two compounds, the surface tensions are grossly underestimated. The estimations by neural networks from σ-moments are not very satisfactory, which suggests that the σ-moments do not contain sufficiently relevant information for our task. In particular, graph machines provide the most accurate estimation as they take into account the structure of the molecule, thereby Difference between measured and predicted ST using graph machines. [START_REF] Andersson | First-Principles Prediction of Liquid/Liquid Interfacial Tension[END_REF] f The root-mean-square error is computed with eq 13 for the 23 molecules of the test set; the estimations for the surface tension are averaged over the 25 models (out of 250) that have the smallest VLOO scores for the multilayer Perceptron and the graph machines (71 and 95 parameters, respectively). g n/a: not applicable.

exempting the model designer from the computation and selection of features.

CONCLUSIONS

Predicting surface tensions of liquids from first-principles is not straightforward due to the complexity of the interactions taking place at gas/liquid interfaces. It remains still inaccessible even by advanced methods such as COSMO-RS able to calculate accurately interfacial tensions between two liquid phases. [START_REF] Andersson | First-Principles Prediction of Liquid/Liquid Interfacial Tension[END_REF] In the present study, the estimation of surface tensions of cosmetic oils was performed by the group contribution method, by the corresponding-states method, by nonlinear regression (multilayer Perceptron) from COSMO-RS σ-moments, and by regression on graphs (graph machines). To the best of our knowledge, this is the first study that compares the accuracy of various methods for estimating the surface tensions of the same molecules and where the models are tested on two sets of molecules; it is also the first report of the use of graph machines for surface tension estimation. The database of surface tensions was created with special emphasis on the diversity of the molecules and on the reliability of the experimental results. It was shown that the group contribution method failed to provide accurate estimations (root-mean-square estimation error = 6.7 mN•m -1 ); the corresponding-states method and nonlinear regression from σ-moments yielded estimations with similar accuracies (root-mean-square estimation errors = 2.4 and 2.7 mN•m -1 respectively); graph machines, which make estimations solely from the molecular structure, provided the most accurate estimations (0.9 mN•m -1 ). Pending the development of methods, based on first-principles, for surface tension estimation only from molecular structures, graph machines are shown to be able to provide, from molecular structures, estimations whose errors are on the order of the experimental uncertainty on the available data. An interactive demonstration of the graph machine computations, based on the Docker free software technology (existing for many computer platforms) is available for download (see the Supporting Information). In the present study, all measurements were taken at the same temperature. As surface tension varies with temperature, we intend to use the ability of graph machines to take into account exogenous variables (e.g., temperature) for prediction. *Phone: 33-3-20336364. E-mail: jean-marie.aubry@univ-lille1.fr (J.-M. Aubry).

ORCID

Francois Duprat: 0000-0002-2889-1701 Veŕonique Nardello-Rataj: 0000-0001-8065-997X

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS

One of the authors (V. Goussard) thanks the company Oleon for financial support by means of doctoral fellowship.

  used an 86-component database containing almost exclusively hydrocarbons (linear and branched alkanes and alkenes, aromatics) to develop the model for the surface tension based on a CSP approach. It includes two reference fluids, methane and n-octane:[START_REF] Zuo | Corresponding-States and Parachor Models for the Calculation of Interfacial Tensions[END_REF]

Figure 1 .

 1 Figure 1. Family distribution (as percentages) for the molecules of the complete (left) and test (right) data sets.

Figure 2

 2 gives an example of the σ-profiles of three cosmetic oils: propylheptyl caprylate, dioctyl ether, and n-tetradecane. None of these compounds has Lewis acidity (hydrogen-bond (HB) donor region), but Lewis basicity (HB acceptor region) is present. The σ-profiles show the difference in the Lewis basicity, propylheptyl caprylate having the strongest one since the curve exhibits the most significant peak above 0.01 e.Å -2 .

Figure 2 .

 2 Figure 2. Plots of σ-surfaces and σ-profiles of three cosmetic oils: propylheptyl caprylate (blue), dioctyl ether (red), and tetradecane (green) using the COSMO-RS model.

Figure 3 .

 3 Figure 3. Variation of the surface tensions (ST) for various families of linear molecules: alkanes (blue circles), ethers (empty purple circles), acetates (red triangles), alcohols (empty orange triangles), and carboxylic acids (yellow crosses) with the number of carbon atoms and with the sum of carbon and silicon atoms for Si-containing compounds (green squares for linear silicones and empty brown squares for cyclic silicones).

a

  RMSTE value (mN•m -1 ) of the model (out of 100) having the smallest RMSTE for the 244 molecules of the training/validation set. b Mean and standard deviation (in parentheses) of the VLOO scores (mN•m -1 ) averaged over the 10 models (out of 100) having the smallest VLOO scores for the 244 molecules of the training/validation set.

Figure 4 .

 4 Figure 4. Surface tension estimation from σ-moments by neural networks. VLOO surface tension (ST) estimates computed by neural network (multilayer Perceptron with 10 hidden neurons) for the 244 molecules of the training/validation set (gray hollow circle) and surface tension estimates for the 25 molecules of the test set (black filled circle) vs measured values of the surface tension. The VLOO estimates for the molecules of the training/validation set and the estimates for the molecules of the test set are computed with the 10 models (out of 100) that have the smallest VLOO scores. The dashed gray line is the regression line for the training/validation set, while the dashed black line is the regression line for the test set.

a

  RMSTE value (mN•m -1 ) of the model (out of 100) having the smallest RMSTE for the 244 molecules of the training/validation set. b Mean and standard deviation (in parentheses) of the VLOO scores (mN•m -1 ) averaged over the 10 models (out of 100) having the smallest VLOO scores for the 244 molecules of the training/validation set.

Figure 5 .

 5 Figure 5. Surface tension estimation from SMILES by graph machines. VLOO surface tension (ST) estimates computed by graph machines (node function with seven hidden neurons) for the 244 molecules of the training/validation set (gray hollow circle) and surface tension estimates for the 25 molecules of the test set (black filled circle) vs measured values of the surface tension. The VLOO estimates for the molecules of the training/validation set and the estimates for the molecules of the test set are computed with the 10 models (out of 100) that have the smallest VLOO scores. The dashed gray line is the regression line for the training/validation set, while the dashed black line is the regression line for the test set.

Table 3 .

 3 Difference between Measured and Predicted Surface Tension (ST) for 23 Cosmetic Oils a Mean measured surface tension values (ST) are reported in mN•m -1 at 23.5 °C. b Difference between measured and predicted ST using group contribution method implemented in the IBSS software. c Difference between measured and predicted ST using corresponding state-based method included in IBSS software. 29 d Difference between measured and predicted ST using neural network model based on the five σ-moments.6e

■

  ASSOCIATED CONTENT * S Supporting Information The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.7b00512. Names, SMILES notations, molecular formulas, CAS RN, σ-moment descriptors, and measured surface tensions values in mN•m -1 at 25 °C of the 244 molecules of the training/validation set, the 25 molecules of the test set, and the 23 molecules of the cosmetic oil set; groupcontribution Conte method estimated, correspondingstates Pitzer's method estimated, corresponding-states Zuo-Stenby method estimated, neural network estimated, neural network VLOO estimated, graph machine estimated and graph machine VLOO estimated values of the surface tension in mN•m -1 for the 269 molecules of the complete set; group-contribution estimated, neural network estimated, and graph machine estimated values of the surface tension in mN•m -1 for the 23 molecules of the cosmetic oil set; Docker Client installation, graph machine demonstration with Docker containers and graph machine results with Docker. (PDF) ■ AUTHOR INFORMATION Corresponding Authors *Phone: 33-1-40794465. Fax: 33-1-40794466. E-mail: arthur. duprat@espci.fr (F. Duprat).

  Kondor et al. used a data set of 1275 surface tension values, collected at different temperatures for 188 molecules, for building and validating the model. When the model was applied to a test set of 225 data points, the squared correlation coefficient R 2 between the estimated surface tensions and the experimental values was equal to 0.96. Stanton et al. described a 10-descriptor model for the prediction of surface tensions of 146 compounds. 8 Kauffmann et al. developed a model based on eight descriptors for a data set of 159 molecules.

•

  Construction of the computational structure: For each acyclic graph, a function is generated by implementing, at each node of the graph, a parametrized nonlinear function, called node f unction (typically an MLP). All node functions are identical and have the same parameters within each graph and for all graphs; therefore, the number of parameters of the resulting model is equal to the number of parameters of the chosen node function. As a result of this construction, the value computed by the output node of each model, which is intended to be an estimate of the property of interest, depends solely on the 2D structure of the molecule and on the values of the parameters of the node function.• Estimation of the parameters of the node function by training from the database: This is done by minimizing the sum of squared errors J(θ) defined in Section 2.2.3. Details of the above steps are provided in a previous paper

Table 1 .

 1 Surface Tension Estimation from σ-Moments by Neural Networks (Multilayer Perceptrons). Estimation of the Quality of Training and of the Predicting Ability for Increasing Neural Network Complexity

	number of parameters	29	36	43	50	57	64	71	78	85	92
	number of hidden	4	5	6	7	8	9	1 0	1 1	1 2	1 3
	neurons
	RMSTE a	2.55	2.41	2.22	1.94	1.75	1.66	1.46	1.28	1.25	1.13
	VLOO score b	3.54 (0.02) 3.29 (0.05) 3.17 (0.09) 2.98 (0.15) 2.96 (0.11) 2.83 (0.21) 2.72 (0.08) 2.71 (0.13) 2.71 (0.12) 2.74 (0.10)

Table 2 .

 2 Surface Tension Estimation from SMILES by Graph Machines. Estimation of the Quality of Training and of the Predicting Ability for Increasing Graph Machine Complexity

	number of parameters	31	44	59	76	95	116	139	164
	number of hidden neurons	3	4	5	6	7	8	9	10
	RMSTE a	1.99	1.33	0.94	0.75	0.49	0.40	0.31	0.18
	VLOO score b	2.77 (0.08)	2.40 (0.15)	2.31 (0.27)	2.05 (0.21)	1.71 (0.29)	1.85 (0.26)	1.90 (0.26)	2.46 (0.77)