Laurent Pagani
email: laurent.pagani@obspm.fr

Evelyne Roueff

Pierre Lesaffre

Ortho-H 2 and the age of interstellar dark clouds

Keywords: ISM, abundances -ISM, clouds -evolution -ISM, molecules -astrochemistry

Interstellar dark clouds are the sites of star formation. Their main component, dihydrogen exists under two states, ortho and para. H 2 is supposed to form in the ortho:para ratio (OPR) of 3:1 and to subsequently decay to almost pure para-H 2 (OPR ≤ 0.001). Only if the H 2 OPR is low enough, will deuteration enrichment, as observed in the cores of these clouds, be efficient. The second condition for strong deuteration enrichment is the local disappearance of CO, which freezes out onto grains in the core formation process. We show that this latter condition does not apply to DCO + , which, therefore, should be present all over the cloud. We find that an OPR ≥ 0.1 is necessary to prevent DCO + large-scale apparition. We conclude that the inevitable decay of ortho-H 2 sets an upper limit of ∼6 million years to the age of starless molecular clouds under usual conditions.

Introduction

The low mass star formation process is relatively well understood, however its first step details remain uncertain, especially regarding the question of time and age. Still much debated questions are: how long it takes for a cloud to form from the diffuse atomic medium, how long a molecular cloud lives, how long it takes to form a condensation that will evolve into a prestellar core that subsequently collapses to form a protostar. Today, much debate exists upon the lifetime of clouds, e.g. Hartmann et al. (2001) defending a short lifetime (1 to a few million years) while Tassis & Mouschovias (2004) and Mouschovias et al. (2006) claim a typical age of 10 million years, but most arguments are either of limited statistical significance or model-dependent (unknown magnetic field strength, small scale turbulence rapid dissipation, etc.). Clues are needed and while the usual chemical age modeling is not satisfying (problem of unknown initial conditions), we present here a simple constraint based on basic chemistry: the absence of DCO + in dark cloud envelopes can only be explained if the clouds are young enough and we propose to determine an upper limit for this age which is as much independent from the initial conditions as possible.

The abundance of CO and DCO + , and the deuterium chemistry

Dark clouds contain dust grains embedded in interstellar gas, itself mostly molecular as these clouds are self-shielded against our Galaxy's UV background. These dark clouds are cold (typically 10 K, with cores as low as 7 K, Pagani et al. 2003Pagani et al. , 2005)), and deep inside, where the density reaches a few 10 4 cm -3 , the heavy species freeze out onto grains to form ices. Both these ices and the gas are subject to a rich chemistry. Recently, strong deuteration (deuterium enhancement with respect to hydrogen carriers: DCO + versus HCO + , HDCO and D 2 CO versus H 2 CO, etc.) has been recognized as a useful tool to better understand both the chemistry itself and the star formation process (Bergin & Tafalla 2007;[START_REF] Ceccarelli | Protostars and Planets V[END_REF]).

H 2 is the main hydrogen reservoir in dark clouds, and similarly, HD is the main deuterium reservoir. Its relative abundance to H 2 is ∼3 × 10 -5 [START_REF] Hébrard | Far Ultraviolet: Five Years of Discovery with FUSE[END_REF]Linsky et al. 2006). Watson (1976) was the first to invoke chemical deuterium enrichment via the important reaction

H + 3 + HD → H 2 D + + H 2 + ∆E (1)
∆E = 232 K when all species (reactants and products) are in their ground state. In cold dark clouds, the forward reaction is strongly favorable, creating H 2 D + /H + 3 ratios much larger than the original HD/H 2 ratio. If CO is abundant (X[CO] ≈ 1.5 × 10 -4 , where X[] denotes the relative abundance of a species with respect to H 2 ), H + 3 reacts 5 times more often with CO than with HD, and forms HCO + , strongly decreasing the abundance of H 2 D + .

Once H 2 D + is formed, it can be enriched in deuterium by reacting with HD again to form D 2 H + and D + 3 as first recognized by Roberts et al. ( 2003) :

H 2 D + + HD → D 2 H + + H 2 + 187K
(2)

D 2 H + + HD → D + 3 + H 2 + 234K (3)
These two reactions are again favorable in the forward direction at low temperatures because of their exothermicity. Then, these H + 3 isotopologues can easily transfer their deuterons to other species:

CO + H 2 D + → 2 3 HCO + + 1 3 DCO + + 1 3 H 2 + 2 3 HD (4) CO + D 2 H + → 1 3 HCO + + 2 3 DCO + + 1 3 HD + 2 3 D 2 (5) CO + D + 3 → DCO + + D 2 (6) N 2 + H 2 D + → 2 3 N 2 H + + 1 3 N 2 D + + 1 3 H 2 + 2 3 HD (7) N 2 + D 2 H + → 1 3 N 2 H + + 2 3 N 2 D + + 1 3 HD + 2 3 D 2 (8) N 2 + D + 3 → N 2 D + + D 2 (9)
When CO is abundant, H 2 D + is much less abundant as noted above and also reacts with CO preferentially, quenching all the other deuteration paths. CO also reacts with species like N 2 H + and N 2 D + to destroy them:

N 2 H + + CO → HCO + + N 2 (10) N 2 D + + CO → DCO + + N 2 (11)
The strong correlation between CO freeze-out onto grains and deuteration (e.g. [START_REF] Crapsi | [END_REF] confirm the chemical models (e.g. Roberts et al. 2003) as sketched above and has led to the conclusion that CO freeze-out was necessary for strong deuteration to occur. There is however one exception to this statement: the deuteration of CO itself can form DCO + . In first approximation, if H + 3 is mainly destroyed by CO, X[H 2 D + ] varies like X[CO] -1 , whereas the production of DCO + is proportional to X[CO] and X[H 2 D + ] (reaction 4). Thus the dependency upon X[CO] cancels out and X[DCO + ] should remain approximately constant. Figure 1 shows the abundance of X[DCO + ], X[H + 3 isotopologues] and X[N 2 D + ] as a function of X[CO] in a steady-state chemical model. The X[H + 3 isotopologues] is the sum of the abundances of all 3 deuterated H + 3 isotopologues, which represent altogether the main deuteration partners of species like CO and N 2 . The model is the same as in Pagani et al. (2009) but CO now runs from an undepleted relative abundance of 1.5 × 10 -4 down to a depletion by a factor 300 (X[CO] = 5 × 10 -7 ), for a cloud density of 10 4 cm -3 at 10 K. When X[CO] starts to drop, X[H + 3 isotopologues] increases faster than proportional (indicated by the -1 slope) because the D 2 H + abundance depends both on the destruction of H + 3 by CO and on the destruction of H 2 D + by CO, and the D + 3 abundance depends also on the destruction of D 2 H + by CO. This explains that at the beginning, X[DCO + ] increases: the deuteration capabilities increase faster than X[CO] decreases. While X[CO] keeps decreasing, the deuteration capabilities reach saturation and X[DCO + ] becomes proportional to X[CO] (traced by the +1 slope). X[N 2 D + ], a representative gas-phase deuterated species, follows the deuteration capabilities of the cloud, as expected. It is clear that DCO + has a behavior different from the other deuterated species, with an abundance varying only by 50% while the CO abundance drops by a factor of 30, and most importantly by reaching its maximum abundance when CO is high and deuteration capabilities are down. Therefore, one would expect DCO + to be present everywhere in the cloud. It is not the case. Indeed, DCO + is, as the other deuterated species, only detected in the depleted prestellar cores on extents similar to NH 3 (Butner et al. 1995;Juvela et al. 2002) and nowhere else. This can only be explained if a high H 2 OPR is present. Indeed the H 2 OPR is an important controller of the cold cloud deuteration chemistry, as already discussed in Pineau des Forêts et al. (1991); Flower et al. (2006) and Pagani et al. (1992Pagani et al. ( , 2009)).

The role of ortho-H 2

H 2 has 2 different possible nuclear spin states due to the proton spin of 1 2 : parallel (ortho, I = 1, weight 2I+1 = 3 ) and anti-parallel (para, I = 0, weight = 1). The ortho state corresponds to the odd rotational levels (J = 1, 3, 5, ), the para state to the even levels. The first ortho state (J = 1) is 170 K above the ground (J = 0) para state. The ortho-para exchange is possible neither via radiative processes, nor via inelastic collisions. It can occur in the gas phase and on the surface of dust grains (Le Bourlot 2000). The conversion rate on a solid surface is however highly uncertain (Watanabe et al. 2010;[START_REF] Ceccarelli | Protostars and Planets V[END_REF]) and it is most probable that this is mainly realized by gas-phase chemical reactions with H + or H + 3 ions (Dalgarno et al. 1973;Flower & Watt 1984;Le Bourlot 1991;Flower et al. 2006). The importance of ortho-H 2 in the deuteration process was first studied by Pineau des Forêts et al. (1991); Pagani et al. (1992), and later on, Flower et al. (2006) revealed its importance as a deuteration regulator. However, they considered the consequences of the OPR evolution only for the prestellar core formation itself. In a study of L183 main prestellar core, Pagani et al. (2009) noted that the H 2 OPR evolution timescale required by the observed N 2 D + /N 2 H + ratio was on the order of the free-fall time. They also noted that the OPR had to drop below 0.01 for N 2 D + to become detectable.

The H 2 OPR value is not easily accessible. H 2 does not emit in cold dark clouds and is seen only in shocked regions, especially protostar outflows, or in absorption in UV in low extinction regions. The energy release from the H 2 formation reaction, assumed to be 1.5 eV in the equipartition hypothesis, is large enough to populate many levels of both ortho and para species: it is hence generally accepted, though not certain, that the OPR is 3:1 when H 2 forms on grains. Once released in the gas phase, H 2 tends to relax towards its thermal equilibrium state but never reaches it because at trace abundance levels, ortho-H 2 destruction is compensated by fresh H 2 formed onto grains from residual H and from destruction of ions such as H + 3 and H 2 D + that can release an ortho-H 2 (Le Bourlot 1991; Flower et al. 2006). It is therefore necessary to model the chemical evolution of the H 2 OPR to understand the evolution of dark clouds. Once H 2 is formed, it can react with H + and H + 3 to exchange states (Le Bourlot 1991;Flower et al. 2006). For example, o-H 2 + H + → p-H 2 + H + + ∆E (12) ∆E = 170 K is the energy released as H 2 goes from J=1 to J=0 levels. The reverse endothermic reaction is thus difficult to obtain in a cloud at 10 K and ortho-H 2 slowly decreases to low abundance levels (≈10 -3 -10 -4 ). Similarly, in reaction 1, except for HD, all the species have ortho and para states and ∆E therefore depends in fact on the considered initial and final states. If only para states are involved in the products, ∆E is maximal (∆E = 232 K with para-H + 3 and 265 K with ortho-H + 3 , HD being always in its ground state, J = 0) but if both products are in their ortho ground state while the reactants are in their lowest state, para-H + 3 ground state and HD (J = 0), ∆E becomes negative and the reaction is slightly endothermic (H 2 D + lowest ortho state is 87 K above the para ground state, and the total required energy is 25 K). Similarly, while the reverse reaction is impossible in a 10 K cloud with both para reactants, it becomes rapid with ortho-H 2 D + and ortho-H 2 even at 10K, as this channel is slightly exothermic (∆E = 25 K) to produce para-H + 3 . A few such channels are opened with ortho-H 2 allowing the efficient destruction of the H + 3 deuterated isotopologues. Ortho-H 2 is thus an important chemical poison to deuteration in dark clouds. Its abundance is therefore critical to control the deuteration enhancement (Flower et al. 2006;Pagani et al. 2009). Because DCO + is not observed everywhere in the clouds outside the depleted cores where it is detected (Butner et al. 1995;Juvela et al. 2002) while its abundance should be maximal, this sets a minimum level for the abundance of ortho-H 2 . In turn, this minimum level gives rise to an upper limit to the age of the cloud, as ortho-H 2 abundance must eventually decrease to low levels as implied by reaction ( 12), and in order to explain the deuteration enhancements seen in the depleted cores, which, in some cold cores, can become huge (amplification up to 10 12 , Lis et al. 2002;van der Tak et al. 2002;Parise et al. 2004).

The chemical model

To address these questions, we used the deuteration network from Roueff et al. (2005Roueff et al. ( , 2007) in which we have included the ortho-para chemistry as described in Pagani et al. (2009), based itself on the work by Walmsley et al. (2004); Flower et al. (2006) and the new rate coefficients of [START_REF] Hébrard | Far Ultraviolet: Five Years of Discovery with FUSE[END_REF]. We have also included analytical approximations of the spin state-dependent dissociative recombination rates of H + 3 isotopologues from the tables in Pagani et al. (2009). The deuteration network includes all simple species (up to 6 atoms) based on C, N, O, and S with up to 5 deuterium substitutions (CD + 5 ). For those reactions with a small endothermicity (E/k B < 600 K), we have corrected the endothermicity when the reaction with ortho-H 2 was involved (e.g., CD + H 2 → CH + HD is endothermic with para-H 2 and exothermic with ortho-H 2 ). For all reactions with rare species producing H 2 (e.g., HCN + + H → CN + + H 2 ), we have always assumed that H 2 is released in its para form, the most favorable energetically (similarly, only ortho-D 2 is considered in similar reactions with deuterated species). This does not noticeably accelerate the decay of ortho-H 2 . We checked that after, e.g., 10 5 years, the para-H 2 production rate is 99.9% dominated by reactions of ortho-H 2 with H + 3 and H + . Our model is not aimed at reproducing accurately a given cloud but it allows to follow the evolution of the H 2 OPR in a molecule rich gas. We have computed the chemical evolution of a cloud of constant density (n(H 2 ) = 10 4 cm -3 ), constant temperature (10 K), average grain radius of 0.1 µm and cosmic ionization rate of 1 × 10 -17 s -1 , which represent the usual conditions met in starless dark clouds. We have started with all species being atomic except H which is considered to be already entirely converted to H 2 and D to HD and we have considered different OPR starting ratios from 3 down to 3 × 10 -3 to account for the possibility of H 2 production OPR below 3. Figure 2 shows the evolution of DCO + and ortho-H 2 with time. We have traced the detectability limit of DCO + for a cloud with a H 2 column density of 10 22 cm -2 , considering that a DCO + J:1→0 line with a width of 0.5 km s -1 and an intensity of 0.1 K would be easily detectable with present day radiotelescopes. With the RADEX non-LTE radiative transfer model (van der Tak et al. 2007), we find that a DCO + column density of ∼2 × 10 11 cm -2 is detectable, equivalent to a relative abundance to H 2 of ∼2 × 10 -11 . Figure 2 shows that for any initial H 2 OPR below 0.1, DCO + should be detected throughout the cloud in less than 3 × 10 5 years. For an initial OPR = 0.1, the DCO + becomes marginally detectable around ∼5 × 10 5 years, disappears and finally becomes clearly detectable after 2 My. For OPR > 0.1, it takes 3 to 6 million years to become detectable, i.e. when the H 2 OPR drops below ∼0.03. Compared to the H 2 OPR drop from 3 to ∼0.01 in depleted cores which takes less than 2 × 10 5 years to happen (Pagani et al. 2009), the difference is twofold: as density increases during the prestellar core formation, 1) the chemistry accelerates, and 2) depletion sets in and H + is no more destroyed by the other species which have now disappeared from the gas phase (H + 3 becomes also more abundant but to a lesser extent). H + can even become the dominant ion. The higher abundance of H + and H + 3 and the accelerated chemistry provoke the rapid decline of ortho-H 2 .

On the same figure, we also trace the abundance of the sum of the deuterated H + 3 isotopologues in the case H 2 OPR = 3. It shows that under 1 My, deuteration is dominated by the carbon chemistry (CH 2 D + , C 2 HD + , etc.) which is not dependent upon the ortho-H 2 abundance due to exothermicities close to 400 K. It is only beyond 1 My that the H + 3 deuteration network takes over when ortho-H 2 starts to be negligible. If ortho-H2 is low since the beginning, then the H + 3 isotopologue contribution to the abundance of DCO + is important from the start which explains the disappearance of the dip between 1 and ∼3 My.

The model is also sensitive to two other initial conditions: the abundance of metals and the cosmic ray ionization rate. Metals, especially alkali and alkaline earth metals like Ca, Na, K, and others like Fe, have a direct influence on the electronic equilibrium of the chemical model which in turn changes the dissociative recombination efficiency for DCO + and therefore its abundance. For standard ζ (1 × 10 -17 s -1 ) and H 2 OPR (3), we have varied their abundance from 3.4 × 10 -8 to 1.3 × 10 -7 . For the lowest abundance (50% of the default one), DCO + becomes detectable much earlier ( 2 × 10 5 y) than in the standard case, drops below the detectability limit after 1 My and is detectable again after 4 My. For the highest abundance (twice the default one), DCO + is detectable after 7.5 My. The cosmic ray ionization rate (Fig. 3): if the rate increases, the production of H + and H + 3 increases, therefore accelerating the decay of ortho-H 2 and increasing the abundance of HCO + and DCO + . For ζ = 3 × 10 -17 s -1 , DCO + becomes detectable in less than 10 5 years. Values of ζ in between 1 and 3 × 10 -17 s -1 would bring our model in agreement with some recent estimates of cloud lifetime (Enoch et al. 2008;Hatchell & Fuller 2008). Conversely, for ζ = 3 × 10 -18 s -1 , the chemistry is considerably slowed down, DCO + becomes detectable only after 28 My and steadystate is reached only after 47 My. It must be noted however that such a low ionization rate seems unlikely following Padovani et al. (2009).

Conclusions

Ortho-H 2 abundance is most probably the main limiting parameter for DCO + production. The absence of DCO + outside depleted cores, while its abundance is expected to be maximum for CO abundances close to 10 -4 , is due to a high ortho-H 2 abundance in dark clouds (OPR ≥ 0.1) maintained long enough. Such high abundance can survive thanks to low H + and H + 3 abundances, remaining low via charge or proton exchanges with many neutral species. However, ortho-H 2 eventually disappears. This sets an upper limit to the age of DCO + -quiet dark clouds of 3-6 million years after H 2 has formed for the most probable initial conditions. Higher cosmic ray ionization rates would increase the abundance of H + and H + 3 and can thus only shorten the maximal possible age of the clouds to values which are too small while lower rates would tend to relax the age constraint but seem very unlikely. This result strengthens the results of Hartmann et al. (2001); Enoch et al. (2008); Hatchell & Fuller (2008) and is not compatible with those of Tassis & Mouschovias (2004) and Mouschovias et al. (2006). From the chemical point of view, cold dark clouds seem to be really short-lived. 3 isotopologues is also traced as a thick gray line (for the case H 2 OPR = 3). The horizontal dotted line in the lower box indicates the limit of detection of DCO + and the circles in the upper box mark the corresponding OPR values. Fig. 3.-Variation of DCO + and ortho-H 2 abundances in a pseudo time-dependent chemical model for different cosmic ray ionization rates from 1 × 10 -18 up to 1 × 10 -16 . The horizontal dotted line indicates the limit of detection of DCO + , the arrows mark the time at which DCO + abundance crosses the detection limit.

Fig. 1 .

 1 Fig. 1.-DCO + , N 2 D + and H +3 isotopologues abundances as a function of CO abundance in a steady state model. Slopes of +1 and -1 are traced to better visualize the variation slope of the different ions. DCO + , on the contrary to other deuterated species like N 2 D + reaches its maximum abundance for high CO abundances.

Fig. 2 .

 2 Fig. 2.-Variation of DCO + and ortho-H 2 abundances in a pseudo time-dependent chemical model for different starting H 2 OPRs from 3 × 10 -3 up to 3. The abundance of the sum of the deuterated H +3 isotopologues is also traced as a thick gray line (for the case H 2 OPR = 3). The horizontal dotted line in the lower box indicates the limit of detection of DCO + and the circles in the upper box mark the corresponding OPR values.

We are grateful to Patrick Hennebelle for fruitful discussions. ER acknowledges the support of the ANR N0. 09-BLAN-020901.


Hoel Le

Capitaine Lina
email: hoel.lecapitaine@univ-nantes.fr

A relevance-based learning model of fuzzy similarity measures

Matching pairs of objects is a fundamental operation in data analysis. However, it requires to define a similarity measure between objects to be matched. The similarity measure may not be adapted to the various properties of each object. Consequently, designing a method to learn a measure of similarity between pairs of objects is an important generic problem in machine learning. In this paper, a general framework of fuzzy logical-based similarity measures based on T-equalities derived from residual implication functions is proposed. Then a model allowing to learn the parametric similarity measures is introduced. This is achieved by an online learning algorithm with an efficient implication-based loss function. Experiments on real datasets show that the learned measures are efficient at a wide range of scales, and achieve better results than existing fuzzy similarity measures. Moreover, the learning algorithm is fast, so that it can be used in real world applications, where computation times are a key-feature when one chooses an inference system.

I. INTRODUCTION

Whilst the similarity is an essential concept in human reasoning, and plays a fundamental role in theories of knowledge, there is not an unique and general-purposed definition of similarity. The reason for this lack of definition comes from the fact that one can find practical cases where similarity properties are not satisfied (e.g. symmetry, indiscernibility or transitivity, see [START_REF] Tversky | Similarity, separability, and the triangle inequality[END_REF]). Indeed, several studies (see [START_REF] Cock | On (un)suitable fuzzy relations to model approximate equality[END_REF], [START_REF] Klawonn | Should fuzzy equality and similarity satisfy transitivity?[END_REF] and references therein) have shown that similarity measures do not necessarily have to be transitive, implying a contradiction with the most usual approach of comparison, based on geometrical assumptions in the feature space.

Fuzzy set theory provides a consistent basis for information processing, and an elegant, mathematically well-founded, representation of the uncertainty in the data. Since the data to be processed are often imprecise, using fuzzy set theory or its derivatives (e.g. possibility theory or belief function theory) has become a common approach in recent years [START_REF] Hullermeier | Fuzzy methods in machine learning and data mining: status and prospects[END_REF].

In this paper, similarity measures are defined by the use of ⊤-equalities derived from fuzzy residual implications. It does not suffer from the drawbacks of the conventional metric approaches, and allows to obtain concave or convex isosimilarity contours.

However, the number of similarity measures induced by the proposed framework is infinite, due to the infinite number of triangular norms (t-norms). In practice, the user must choose the similarity measure, and this wide range of choice is problematic. Consequently, learning the similarity measure from the available data is proposed. To this aim, the similarity measure is defined such that relevance degrees between objects are respected when then are ranked according to their pairwise similarity.

Additionally, an online learning setting is adopted, where, in contrast to batch methods, samples are considered one at a time. This enable to treat large data sets while keeping satisfying performances [START_REF] Chechik | Large scale online learning of image similarity through ranking[END_REF].

This paper is organized as follows. Section II first reviews the main approaches dealing with similarity measurement. Then, ⊤-equalities are used to define new classes of (parametric) similarity measures. Section III presents the learning algorithm and its properties, as well as examples. In Section IV, the use of the learning algorithm for supervised classification is described, and some comments on performances are given. Finally, conclusion and perspectives are drawn in Section V.

II. FUZZY SIMILARITY MEASURES

A. Basic material

Aggregating values plays an important role in decisionmaking systems. Given n values, an aggregation operator is a mapping A : [0, 1] n → [0, 1] satisfying boundary conditions and monotonicity. In the literature, one finds many aggregation operators, e.g.: t-norms , OWA (Ordered Weighted Averaging) operators, γ-operators, or fuzzy integrals. They belong to several categories, depending on the way the values are aggregated: conjunctives, disjunctives, averaging, and mixed operators. The interested reader can refer to [START_REF] Calvo | Aggregation operators: Properties, classes and construction methods[END_REF], [START_REF] Grabisch | Aggregation Functions, ser. Encyclopedia of Mathematics and its Applications[END_REF] for large, yet comprehensive, surveys on aggregation operators.

A t-norm is an increasing, associative and commutative mapping ⊤ : [0, 1] 2 → [0, 1] satisfying the boundary condition ⊤(x, 1) = x for all x ∈ [0, 1]. The most popular continuous t-norms are the minimum ⊤ M (x, y) = min(x, y), the product ⊤ P (x, y) = x y and the Łukasiewicz t-norm ⊤ L (x, y) = max(x + y -1, 0). Various parametric families involving a real value λ lying in a specified domain have been introduced. The parametric t-norms that are used in the sequel are given in Table I.

(x, y) =      Drastic t-norm if λ = ∞ 0 if λ = x = y
x y λ+(1-λ) (x+y-x y) if λ ∈ [0, ∞[ and (λ, x, y) = (0, 0, 0)

Dombi ⊤D(x, y) =          Drastic t-norm if λ = 0 ⊤M (x, y) if λ = ∞ 1 + 1-x x λ + 1-y y λ 1/λ -1 if λ ∈]0, ∞[ Yager ⊤Y (x, y) =      Drastic t-norm if λ = 0 ⊤M (x, y) if λ = ∞ max 1 -(1 -x) λ + (1 -y) λ 1/λ , 0 if λ ∈]0, ∞[ Frank ⊤F (x, y) =          ⊤M (x, y) if λ = 0 ⊤P (x, y) if λ = 1 ⊤L(x, y) if λ = ∞ log λ 1 + (λ x -1)(λ y -1) λ-1 if λ ∈]0, 1[∪]1, ∞[ Schweizer- Sklar ⊤SS(x, y) =              ⊤M (x, y) if λ = -∞ ⊤P (x, y) if λ = 0 Drastic t-norm if λ = ∞ max x λ + y λ -1, 0 1 λ if λ ∈] -∞, 0[∪]0, ∞[
A general problem in fuzzy logic is to handle conditional statements if x, then y where x and y are fuzzy predicates. A widely used method consists in managing them by using functions I : [0, 1]×[0, 1] → [0, 1] such that the truth value of I depends on the initial propositions x and y. We generally speak about an implication function if I is non-increasing in the first variable, non-decreasing in the second variable, I(0, 0) = I(1, 1) = 1, and I(1, 0) = 0, see [START_REF] Mas | A survey on fuzzy implication functions[END_REF] and [START_REF] Baczynski | Fuzzy Implications, ser. Studies in Fuzziness and Soft Computing[END_REF] for recent surveys on fuzzy implication functions.

In this paper, residual implication functions are considered. Given a t-norm ⊤, its corresponding residuum is defined by

I ⊤ (x, y) = sup t {t ∈ [0, 1]|⊤(x, t) ≤ y}. (1)
In the sequel, the following notation is adopted:

• X = {x 1 , • • • , x n } is the (supposed finite) universe of discourse, • C(X)
and F(X) are the sets of all crisps and fuzzy sets in X, respectively, • f A (x), ∀x ∈ X, is the membership function of a fuzzy set A over X.

There are several ways to compare fuzzy values or fuzzy quantities. The first one is based on a broad class of measures of equality based on a distance measure which is specified for membership functions of fuzzy sets. The second category involves set-theoretic operations for fuzzy sets (fuzzy intersection, union, cardinality) [START_REF] Tversky | Features of similarity[END_REF], [START_REF] Bouchon-Meunier | Towards general measures of comparison of objects[END_REF]. Finally, a third way of defining a similarity measure consists in using logical concepts of fuzzy implication, as first suggested in [START_REF] Hirota | Matching fuzzy quantities[END_REF], following the seminal paper of [START_REF] Bandler | Fuzzy power sets and fuzzy implication operators[END_REF]. Note that a fourth approach, relying on morphological operators, has been proposed in [START_REF] Bloch | On fuzzy distances and their use in image processing[END_REF].

First of all, the basic definition of a fuzzy similarity measure is recalled.

Definition 1. A mapping S : F(X) × F(X) → [0, 1] is called a similarity measure if it satisfies (P1) S(A, B) = S(B, A), ∀A, B ∈ F(X). (P2) S(A, A) = 1, ∀A ∈ F(X). (P3) S(D, D c ) = 0, ∀D ∈ C(X). (P4) ∀A, B, C ∈ F(X), if A ⊆ B ⊆ C, then S(A, C) ≤ S(A, B) ∧ S(B, C) or, equivalently ∀A, B, C, D ∈ F(X), if A ⊆ B ⊆ C ⊆ D, then S(A, D) ≤ S(B, C)
However, most of these properties, if not all, are subject to criticisms and debates. Therefore, it contributes to emphasize the lack of a clear definition of a general-purposed similarity measure. The symmetry property (P1) is still subject to experimental investigations: if S(x, y) is the answer to the question how is x similar to y?, then one focus more on the feature x than on y. This corresponds to the notion of saliency [START_REF] Tversky | Features of similarity[END_REF] of x and y: if y is more salient than x, then x is more similar to y than vice versa, which is experimentally confirmed by observing asymmetries in confusion matrices. The property (P2) is also discussed, since in some cases, the similarity of a point to itself is not 1. The property (P4) is the most debatable one. In particular Tversky and Gati [START_REF] Tversky | Similarity, separability, and the triangle inequality[END_REF] present and criticize the segmental additivity, i.e. d(A, C) = d(A, B) + d(B, C) property of metrics. This assumption is rather intuitive, and include a wide class of distance functions: all the Minkowski metrics, and Riemannian curved geometries. When dealing with fuzzy sets, and considering ⊤-equivalences and ⊤-equalities, the triangle inequality is replaced by the ⊤-transitivity.

The metric approach takes its roots from studies on how to measure the distance between two real functions. The basic concept is to consider a fuzzy set as a point in a vector space. The general form of a Minkowski r-metric defined for r ≥ 1 is usually taken:

d r (A, B) = x∈X |f A (x) -f B (x)| r 1/r (2)
This metric induces well known distance functions; Hamming (or Manhattan) for r = 1, Euclidean for r = 2, Tchebychev (or sup distance) for r = ∞ (which can be written as sup x∈X |f A (x) -f B (x)|). Note that for r < 1, d r does not define a metric, since it violates the triangle inequality. The Hamming and Euclidean distances are also denoted L 1 and L 2 norms, respectively. Increasing the value r gives more weight to large differences in feature values.

Starting from a distance function, several methods have been proposed to obtain a similarity measure. The most natural way is to use a non increasing function g such that S(A, B) = g(d(A, B)). This idea relies on the link between a distance in a vector space and a similarity. Intuitively, when the distance increases, the similarity should decrease. Popular choices of this function g are the 'Cauchy-like' function g(x) = 1/(1 + x) suggested by Zimmerman in [START_REF] Zimmermann | Latent connectives in human decision making[END_REF]. Another function, coming from the work of Shepard [START_REF] Shepard | Toward a universal law of generalization for psychological science[END_REF], consists in taking g(x) = exp(-αx). The use of the exponential law comes from the underlying observation that there is a non linear relationship between the distance and the similarity which is concave upward. However, obtaining similarity measures with the help of distance measures implies that the similarity satisfies the metric properties, which is far from obvious in practice. This has led researchers to prefer non-metric models of similarity.

The second way to compare fuzzy values comes from some basic set-theoretic considerations where union, intersection and complementation are defined for fuzzy sets. While metric based measures can be interpreted as proximity of fuzzy sets, the set-theoretic measures can be viewed as an approximate equality. Probably the most famous set-theoretic measure is the consistency index, defined by the supremum over X of A ∩ B. The interested reader can refer to [START_REF] Cross | Similarity and compatibility in fuzzy set theory: assessment and applications[END_REF] for more details on set-theoretic comparison measures.

The last main approach consists in considering the implication degrees of the elements belonging to A over the elements of B [START_REF] Bandler | Fuzzy power sets and fuzzy implication operators[END_REF], and vice-versa. The implication degree is obtained by using one of the fuzzy implication functions described in the previous section. Formally, having x an element of A, and y an element of B, their implication degree is given by I(x, y). More details on related logical comparison measures are given in [START_REF] Le Capitaine | Towards a unified logical framework of fuzzy implications to compare fuzzy sets[END_REF]. In order to obtain the implication degrees of A over B and B over A, we use the bi-implication bI, defined by bI(x, y) = min(I(x, y), I(y, x))

In [START_REF] Baets | Metrics and t-equalities[END_REF], the authors showed that in the case of residual implications, bI is a ⊤-equality if and only if ⊤ is a leftcontinuous t-norm. In this paper, residual implications are used, since the others (S, QL, D) do not define ⊤-equalities.

Theorem 1. Let bI ⊤ be a bi-residual implication function defining a ⊤-equality E ⊤ . For arbitrary A, B ∈ F(X), let

S(A, B) = n A i=1 E ⊤ (f A (x i ), f B (x i )) (4)
for all x i in X, where A is an aggregation operator. Then S is a similarity measure.

Proof: (P1) By definition, E ⊤ (x, x) = 1 holds, for any x ∈ [0, 1]. By boundary conditions on A, the equality

S(A, A) = 1 is obtained. (P2) by commutativity of ⊤-equality, S(A, B) = n A i=1 E ⊤ (f A (x i ), f B (x i )) = n A i=1 E ⊤ (f B (x i ), f A (x i )) = S(B, A) (P3) by definition, I ⊤ (1, 0) = 0, so that E ⊤ (1, 0) = 0. By boundary conditions on A, S(D, D c ) = 0. (P4) since A ⊆ B ⊆ C ⊆ D, for all x i ∈ X, f D (x i ) ≥ f C (x i ) (5) f B (x i ) ≥ f A (x i ) (6)
hold. By non-increasingness in the first variable and non-decreasingness in the second variable of I ⊤ , for all x i ∈ X, we have [START_REF] Calvo | Aggregation operators: Properties, classes and construction methods[END_REF]. Using E ⊤ (x, y) = I(max(x, y), min(x, y)) and (5-6), we obtain

I ⊤ (f D (x i ), f A (x i )) ≤ I ⊤ (f C (x i ), f A (x i )) by Eq. (5) and I ⊤ (f C (x i ), f A (x i )) ≤ I ⊤ (f C (x i ), f B (x i )) by Eq.
E ⊤ (f D (x i ), f A (x i )) ≤ E ⊤ (f C (x i ), f A (x i )) and E ⊤ (f C (x i ), f A (x i )) ≤ E ⊤ (f C (x i ), f B (x i )). Last, monotonicity of A gives S(A, D) ≤ S(B, C) which concludes the proof.
Remark 1. The minimum operator is used in [START_REF] Klawonn | Should fuzzy equality and similarity satisfy transitivity?[END_REF], but any tnorm also fulfill the desired properties, because for any I ⊤ , x ≤ y ⇒ I ⊤ (x, y) = 1 (by ordering property, see [START_REF] Mas | A survey on fuzzy implication functions[END_REF]).

B. Examples

Taking particular t-norms and aggregation operators enables to retrieve well-known similarity measures. For instance, taking the arithmetic mean and ⊤ P , ⊤ L , the two measures proposed in [START_REF] Wang | New similarity measures on fuzzy sets and on elements[END_REF] are retrieved. Remark that for the arithmetic mean and the Łukasiewicz implication, the exact contraposition of the L 1 norm is obtained, and taking the minimum instead of the mean gives the contraposition of the Tchebychev norm. More details on fuzzy similarity measures that can be obtained with the generic framework are given in [START_REF] Le Capitaine | Towards a unified logical framework of fuzzy implications to compare fuzzy sets[END_REF].

Remark 2. In this paper, the arithmetic mean, the minimum and maximum operators are taken as overall aggregation operators . However, any aggregation operator such as geometric means, OWA operators or fuzzy integrals can be used. The use of the latter is currently under study in order to take into account feature dependencies, as in [START_REF] Santini | Similarity measures[END_REF].

For illustration purpose, Figure 1 shows some examples of fuzzy similarity measures as well as the influence of the parameter λ for the Hamacher and the Dombi-based measures. The different plots show the similarity value of a given fuzzy set A = {0.5/x 1 , 0.5/x 2 } to all the possible two-dimensional fuzzy sets B for various I ⊤ . As one could expect, the closer to x 1 or x 2 , the higher the similarity. One can also note that different I ⊤ lead to different shapes of iso-similarity contours. According to the third line of Figure 1, one can remark that the Hamacher similarity measure gives a concave contour with a low λ, which becomes convex when the parameter increases. Consequently, the similarity measure is much more flexible than usual metrics, since convex and concave shapes can be modelized. The interesting point is that a measure without segmental additivity can be obtained, which is consistent with the axioms of Tversky and Gati, by choosing a parameter making the measure concave. Another parameter value allows to obtain convex iso-similarity contours, yielding a usual metric satisfying the triangle inequality. This property is very important, since a concave contour of iso-similarities means that the triangle inequality is violated [START_REF] Shepard | Attention and the metric structure of the stimulus space[END_REF]. However, when dealing with fuzzy sets, the ⊤-transitivity is preserved, as discussed earlier.

Moreover, the measures enable to capture the idea that similarity is easier to quantify and makes more sense locally (i.e. small variations) than far away in the feature space, where comparisons and judgments of similarity are difficult. An appealing property of parametric similarity measures is that it allows to rank objects in a different order depending on the parameter λ. Keeping the Hamacher and Dombi based similarity measures, three fuzzy sets A, B and C are used in this experiment. They are defined by:

• A = {0.7/x 1 , 0.05/x 2 , 0.32/x 3 , 0.07/x 4 , 0.10/x 5 }, • B = {0.82/x 1 , 0.75/x 2 , 0.36/x 3 , 0.90/x 4 , 0.04/x 5 }, • C = {0.45/x 1 , 0.34/x 2 , 0.69/x 3 , 0.57/x 4 , 0.16/x 5 }. They can be visually inspected in Figure 2. Even for a human, the ranking of S(A, B) and S(A, C) is not an easy task. In Table II, the amount of similarity between A, B and C with various λ values are given. In this table, maximum similarity values (with respect to the other one) are in bold font. According to this table, it can be seen that one can find a value of λ 1 such that S(A, B) > S(A, C) and a value λ 2 such that S(A, B) < S(A, C). Consequently, the three objects A, B and C are ranked differently depending on λ. If one cannot find two values of λ such that the ranking of three objects is different, then no matter how the measures are learned, they are all equivalent in terms of information retrieval, and do not provide efficient similarity measures.

III. LEARNING THE SIMILARITY FUNCTION

A. Loss function

In this paper, an online learning procedure is adopted. In this setting, the algorithm sequentially receives samples, and predicts an output. Once the output is obtained, the algorithm gets a feedback indicating its goodness. Afterwards, parameters can be changed so that the probability of correct output increases in the next step. An appealing property of online algorithms is that they are relatively simple to implement and quickly (i.e. with a few number of iterations) provides good performances [START_REF] Chechik | Large scale online learning of image similarity through ranking[END_REF]. Moreover, the learning algorithm does not require to use all the learning set, so that a great improvement in terms of computation is obtained when one considers large datasets.

Given information related to the relevance degree of matching between objects described by the help of fuzzy sets, the problem of learning a similarity function S is addressed. A relevance degree R can be understood as a pairwise function of objects a and b which states how strong a and b are related. For instance, in supervised classification, this could encode that both are in the same class. The relevance degree may be defined by user's input, which may be based on results of queries, or by knowledge on objects (supervised classification). More formally, there is a set C of n objects and some relevance degrees R between objects belonging to C. The relevance degrees take their values in the unit interval, 1 meaning complete satisfaction, whereas 0 means that there is no reason at all to match the two objects. Supposing that R(a, b) > R(a, c), the aim is to define and learn a parametric similarity measure S λ such that

S λ (a, b) → S λ (a, c) < 1 (7)
The → operator is a fuzzy implication. In the sequel, only residual implications are considered but this can be adapted to any fuzzy implication. The hinge loss function for every triplet a, b, c is defined by:

ℓ λ = max 0, S λ (a, b) → S λ (a, c) -1 (8)
Naturally, the proposed learning framework can be built with many other fuzzy implication functions (parametric included). For simplicity, the Łukasiewicz implication is used in the sequel, so that the loss becomes

ℓ λ = max 0, -S λ (a, b) + S λ (a, c) (9)
The aim is to minimize the total loss L λ on the learning set. The total loss is defined by summing up the individual losses over all the triplet of the learning set:

L λ = (a,b,c)∈C ℓ λ (a, b, c) (10)
However, even in the case of moderately large datasets, the number of all possible triplets is very large and exhaustive computation becomes intractable in practice. Therefore, an online learning scheme is adopted. Consequently, the minimum loss is not searched in the entire learning set, but with randomly selected samples of the set. More precisely, the sample a is randomly selected in the learning set, and b, c are also uniformly sampled such that a and b share the same class, while c belongs to another class.

B. Optimal updating

Since the similarity measure is determined by its parameter value λ, its optimal value is searched by using an online learning algorithm based on sequential updates of λ. Depending on the similarity measure, the initial value λ 0 varies, see Section IV for comments, and subsection II for a discussion on the various similarity measures that are used in the sequel. It leads to find λ such that:

λ i = argmin λ 1 2 λ -λ i-1 2 + αℓ λ (a, b, c) , (11)
where α ≥ 0. In other terms, during the optimization process, λ i is selected in order to obtain a trade-off between minimizing the loss on (a, b, c) and staying quite close to its previous value λ i-1 . This trade-off is controlled by the 'aggressiveness' parameter α. Naturally, if α is set to zero, then the optimal λ i value is equal to λ i-1 . In contrast, setting α to a high value imposes an important weight on the loss function. It is clear that when the loss ℓ λ (a, b, c) is equal to zero, then the optimal parameter value does not change from the previous iteration, i.e. λ i = λ i-1 . Otherwise, the objective function L is defined by

L(λ) = 1 2 λ -λ i-1 2 + α -S λ (a, b) + S λ (a, c) (12)
The optimal solution with respect to λ is such that the gradient of L, ∂L(λ)/∂λ, vanishes:

∂L(λ) ∂λ = λ -λ i-1 -α ∂ S λ (a, b) -S λ (a, c) ∂λ = 0 (13)
Let G i be the gradient value

∂ S λ (a,b)-S λ (a,c) ∂λ
. Therefore, the optimal new value λ is given by

λ = λ i-1 + αG i
The corresponding algorithm is described in Algorithm 1.

Naturally, the value G i depends on the similarity measure that is used. The corresponding gradients for each similarity measure are given in Table III. For writing convenience,

f 1 = max(a k , b k ) and f 2 = min(a k , b k ) for k ∈ N .
The performance of the learning algorithm also depends on the initial value λ 0 , which is discussed in Section IV-C. When λ is out of definition bounds of the t-norm (e.g. negative values for Hamacher one), the last convenient value λ i-1 is returned. An important case is the one of bound values. When λ reaches these bounds, then the equivalence obtained with the corresponding t-norm is used. For instance, if λ = 0 for the SS norm, then the equivalence obtained by the product is used, thanks to the continuity around 0 of this t-norm. Note also that for a reasonable number of iterations (e.g. n 2 /2), limit bounds such as infinity are not reached. For illustration purpose, three

S ∂ ∂λ S λ (a k , b k ) S H (f 2 -f 1 f 2 )(f 1 -f 2 ) f 2 (λ i + f 1 -λ i f 1 ) + f 1 -f 2 2 S Y -(1 -f 2 ) λ -(1 -f 1 ) λ 1/λ 1 λ 2 λ log(1 -f 2 )(1 -f 2 ) λ -log(1 -f 1 )(1 -f 1 ) λ (1 -f 2 ) λ -(1 -f 1 ) λ -log((1 -f 2 ) λ -(1 -f 1 ) λ ) S D 1 -f 2 f 2 λ - 1 -f 1 f 1 λ 1/λ     - 1 λ 2 log 1 -f 2 f 2 λ - 1 -f 1 f 1 + 1 λ log 1-f 2 f 2 1-f 2 f 2 λ -log 1-f 1 f 1 1-f 1 f 1 λ log 1-f 2 f 2 λ -1-f 1 f 1     S F 1 log(λ) 2    log(λ) 1 + (λ f 2 -1)(λ-1) (λ f 1 -1) (λ f 1 -1)(λ f 2 + f 2 λ f 2 -λ f 2 (f 2 /λ) -1) -(λ f 2 -1)(λ -1)λ f 1 f 1 (λ f 1 -1) 2 - 1 + (λ f 2 -1)(λ-1) (λ f 1 -1) λ    S SS 1 + f λ 2 -f λ 1 1/λ - 1 λ 2 log(1 + f λ 2 -f λ 1 ) + 1 λ log(f 2 )f λ 2 -log(f 1 )f λ 1 1 + f λ 2 -f λ 1 1 1.5 2
G i ← ∂ S λ (a,b)-S λ (a,c) ∂λ ⊲ See Table III. 6: λ ← λ i-1 + αG i 7:
until convergence 8:

return λ 9: end procedure fuzzy sets are used for which the ranking of the similarities S(A, B) and S(A, C) is natural. Such fuzzy sets are plotted in Figure 3. Here, there are only three samples, so they are repeatedly selected into the learning algorithm, where α = 1. If the difference between two similarities increases, then the similarity measure is more efficient for the discrimination of other objects belonging to the learning data (and hopefully on the test set). The difference S(A, B) -S(A, C) is studied as a function of the iterations. The corresponding graph is plotted in Figure 4 (dashed line). Naturally, since B is more similar to A than C, even at the beginning of the algorithm, the difference is quite large. The most interesting point is the evolution of the parameter λ that makes A and B more and more similar, and A and C more and more dissimilar as λ is updated. This difference does not exceed an upper bound reached around the 130-th update of λ. Now, the three fuzzy sets introduced in the previous section are considered. As already mentioned, the ranking of S(A, B) and S(A, C) is not easy. Consequently, the learning algorithm is run a first time, with the supposition that B is more similar to A than C. The corresponding difference curve is plotted in Figure 4 as a solid line. Then, the algorithm is run a second time with the supposition that C is more similar to A than B. The S(A, C) -S(A, B) difference curve is plotted in Figure 4 as a dotted line. As can be seen, supposing that B is more similar to A than C leads to quickly reach the upper bound (around the 200-th iteration). In contrast, supposing that C is more similar to A than B requires much more iterations for the learning algorithm to reach its upper bound. Note that at the beginning, the difference is negative, that is to say S(A, B) > S(A, C), and as λ is updated, S(A, B) < S(A, C) is obtained, accordingly to the assumption that C is more similar to A than B. This example shows the ability of the learning algorithm to adapt the similarity measure to the data, since the ranking of the similarity between three objects can be modified according to prior knowledge.

IV. EXPERIMENTS

In this section, an application in supervised classification of the proposed learning model is presented.

A. Protocol

Let X = {(x 1 , y 1 ), • • • , (x n , y n )} be a learning set having c classes. Each sample x i is represented by a set of p features, and its discrete class label y i ∈ {1, • • • , c} is known. In order to use the similarity measures previously introduced, each sample x must be described by a fuzzy set. In the sequel, an automatic fuzzification scheme is proposed. Each sample is described by a discrete fuzzy set of c * p elements. The basic idea of the fuzzification is to estimate, for each feature of each class, its mean and standard deviation, giving two matrices M and S of size (c × p). For each sample, the fuzzy set is obtained by evaluating a membership function on each of its features with respect to all classes. For simplicity, the Gaussian membership function is used. It is defined by

f (x|σ, m) = exp - (x -m) 2 2σ 2 ( 14
)
where m is the mean and σ the standard deviation. The general procedure is described in Algorithm 2, where the set X i j is composed of the j-th feature of the class i.

Example 1. For illustration purpose, we give an example for a dataset composed of a mixture of three normal 1-dimensional distributions of 100 points each. The means are equal to -2, 1, 3 and standard deviations are equal to 0.4, 0.9 and 0.7, respectively. Each sample is then described by three membership degrees corresponding to the three classes. The membership functions obtained are plotted in Figure 5.

During the learning phase, 3 samples are randomly selected in each iteration step . The first one, a, is the reference one. The sample b is randomly chosen so that it belongs to the same class of a, and c is randomly selected such that its class is not the class of a. Since a and b are in the same class, and c belongs to another class, R(a, b) > R(a, c) holds. The performance of all similarity measures is evaluated by using standard ranking precision measures based on nearest neighbors. For each sample, all other test samples are ranked according to their similarity to the sample. The number Algorithm 2 Fuzzification of each sample x 1: procedure FUZZIFICATION(X) 2:

Set p ← number of features 3:

Set c ← number of classes 4:

for i = 1 to c do 5:

for j = 1 to p do 6:

M (i, j) ← componentwise average(X i j )

7:

S(i, j) ← componentwise standard dev.(X i j )

8:

end for 9:

end for 10:

for each sample x k ∈ X do 11:

for i = 1 to c do 12:

for j = 1 to p do 13:

Set membership degrees of x k using ( 14), 14: with M (i, j) and S(i, j).

15:

end for 16:

end for 17:

end for 18: end procedure of same-class samples among the top k similar samples is computed, giving the precision measure for this sample. When averaged over all samples, an average precision (AP) measure as a function of k is obtained:

AP (k) = 1 n n i=1 1 k k j=1 (y(x i ) = y(x j ))
where x j , j = 1, • • • , k are the k most similar objects with respect to x i . The mean average precision (mAP) is obtained for each similarity measure by averaging the precision at level k across all k values:

mAP = 1 k max kmax k=1 AP (k)
Both AP and mAP are commonly used evaluation measures of similarity based classifiers [START_REF] Chen | Similarity-based classification: Concepts and algorithms[END_REF], [START_REF] Chechik | Large scale online learning of image similarity through ranking[END_REF].

B. Datasets

To validate the efficiency of the proposed learning model, a comparison of its performance on various real datasets available at http://archive.ics.uci.edu/ml/ is conducted. The following datasets are considered:

• the well-known IRIS dataset [START_REF] Anderson | The irises of the gaspe peninsula[END_REF], [START_REF] Bezdek | Will the iris data please stand up?[END_REF] The two classes correspond to radar returns showing evidence of some type of structure in the ionosphere, and those that do not. • The Tae dataset consists of evaluations of 151 teaching assistant. Each assistant is described by 5 attributes. The teacher may belong to three classes: low, medium or high. • The Ecoli dataset contains 336 observations described by 7 attributes. Each sample is a sequence of protein, and the 8 classes correspond to their respective localization sites.

• the Liver Disorders dataset contains 345 samples. Each sample constitutes the record of a single male individual, and is described by 7 features such as blood tests, or number of drinks per day. The class, 0 or 1, denotes some sort of liver disorders. • the Newthyroid dataset is composed of 215 instances. The instances are described by 5 chemical features designed to the diagnostic of the functioning of the thyroid gland. This functioning may be normal, hypo or hyper, and corresponds to the three classes of the problem. • the Vowel dataset is composed of 528 samples described by 10 features. The goal is to recognize one of the eleven vowels of British English. • the original Breast Cancer Wisconsin dataset contains 699 samples, described by 10 attributes. Each sample has one of two possible classes: benign or malignant.

As can be seen, the datasets that are considered present a large variety of problems, ranging from linearly separable (e.g. Wine) to hard classification problems (e.g. Yeast).

C. Results

In this section, the learning algorithm of fuzzy similarity measure is compared to standard fuzzy similarity measures issued from the literature. Here, the aim is to analyze the improvement (if any) of the learning scheme compared to usual similarity measures. In all experiments, and for all similarity measures, the initial value of λ is set to 10. The influence of this initialization is discussed in the sequel. First, a detailed study on the IRIS dataset is provided. The number of top k samples varies from 1 to k max = 10, and the individual average precision measures for each similarity measure are obtained. The corresponding performances are reported in Table IV, where the best score for each k is reported in bold font. The parametric measures are shown to be the best measures over all top k samples. The Hamacher measure performs particularly well, ranking first in 7 out of 10 levels. The second experiment demonstrates the efficiency of the approach on the twelve aforementioned datasets. Here again, the average precision measure over the top k samples from 1 to 10 is considered, and the mean average precision is obtained by taking the means across the 10 average precision measures. The results for each dataset and each similarity measure are given in Table V, where best scores for each dataset are reported in bold font. The last column is the average rank of each similarity measure over all datasets. For comparison purpose, results obtained with a recent online similarity learning algorithm (OASIS, see [START_REF] Chechik | Large scale online learning of image similarity through ranking[END_REF]) are also given in the last row of the table. In the experiments, the OASIS algorithm uses the same feature vectors as the other similarity measures. Here, the aim is to compare fuzzy similarity measures, so that OASIS is not taken into account for the rank computation. According to Table V, the following remarks can be made.

• Whatever the datasets, S M leads to the worst rank. This measure is a point-wise measure in the sense that it uses a single degree of membership to determine their value. • Although S mL is also a point-wise measure, the average performance of this measure is better than S M . The reason is that the measure uses a combination of two degrees of membership, instead of one for S M . • Whatever the datasets, there are at least two parametric similarity measures that perform better than commonly used similarity measures. • The average rank of each similarity measures gives an overview of their performance compared to the others. The parametric similarity measures can be ranked as follows: S SS ≻ S H ≻ S Y ≻ S D ≻ S F . It is not surprising to observe that the Schweizer-Sklar based measure performs the best results. This is the only t-norm (among the considered ones) that is equal to the four basic t-norms, depending on the λ value. Consequently, the similarity measure derived from this t-norm is more flexible than the others. • When compared to OASIS, parametric equivalences behaves favorably. In particular, all similarity measures except Dombi and Frank give better results than OASIS.

In order to compare multiple similarity measures over multiple datasets, a combination of a Friedman test and a Nemenyi post-hoc test is used, following the recommendations of [START_REF] Demsar | Statistical comparisons of classifiers over multiple data sets[END_REF].

Let R i j be the rank of the j-th similarity measure on the ith dataset. The Friedman test compares the average ranks R j over all datasets (last column of Table V). Under the nullhypothesis, stating that two similarity measures are equivalent, their ranks should be equal (here R j = 5 for all j). The Friedman statistic is given by

χ 2 F = 12N ns(ns + 1)   j R 2 j - ns(ns + 1) 2 4   (15)
where N , the number of datasets, and ns the number of similarity measures are big enough, typically N > 10 and ns > 5. A derived and better statistic proposed in [START_REF] Iman | Approximations of the critical region of the friedman statistic[END_REF] is given by

F F = (N -1)χ 2 F N (ns -1) -χ 2 F ( 16
)
The Friedman test χ 2 F = 68.95 proves that the average ranks are significantly different from the mean Rank R j = 5 expected under the null hypothesis. Moreover F F = 28.04 is distributed according to the F distribution with 9 -1 = 8 and (9 -1) × (12 -1) = 88 degrees of freedom. The pvalue computed by using the F (8, 88) distribution is almost zero, so that the null hypothesis is rejected at a high level of confidence.

If the null hypothesis is rejected, the Nemenyi post-hoc test is proceeded. The performance of two similarity measures is significantly different if the corresponding average ranks differ by at least the critical difference, defined by

CD = q α nc(nc + 1) 6N , (17)
where q α values are based on the Studentized range statistic divided by √ 2, (see [START_REF] Demsar | Statistical comparisons of classifiers over multiple data sets[END_REF] for details). Finally, a (#similarity measures × #similarity measures) matrix that summarizes the results is obtained. Each entry of the matrix is 1 if the difference of ranks is significant, and 0 otherwise. In order to provide a more informative visualization, a new matrix is created where each entry ({i}, {j}) shows the difference of individual ranks R({i}) -R({j}) obtained with similarity measures {i} and {j}, or black if the difference is not statistically significant under the Nemenyi test, (see Figure 6).

According to Figure 6, the following remarks are drawn. Two main sets of measures : {1, 2, 3, 4} and {5, 6, 7, 8, 9} can be identified. They correspond to the commonly used similarity measures and the new parametric similarity measures, respectively. Within the first set, one must distinguish the measure {2}, i.e. S L , which is not different from the second set at significance level α = 0.05. However, at α = 0.10, it becomes significantly worse than the measure S SS . In the second set, two measures can be distinguished. The first is {8}, i.e. S F . In terms of statistical significance, this is the worst parametric measure. Although its average rank is better than those of {2, 4}, this is not significant at level 0.05. One may argue that at significance level α = 0.05, the post-hoc test is not powerful enough to detect any significant differences between the measures. The second interesting measure is {9}, i.e. S SS , which is significantly better than the entire first set. Considering α = 0.10, one may conclude:

• S SS is significantly better than all commonly used similarity measures, • none of the parametric measures is significantly better than the others, • S P , S M and S mL are significantly worst than the parametric measures, • we cannot conclude on the statistical significance of S L , except for S SS (worse) and S M (better). The next experiment is the analysis of the initial value of λ with respect to the performance obtained with the learning algorithm, where α = 1. In order to make the reading clear, only the results for the Hamacher based similarity measure are reported, although similar comments are valid for the other parametric measures. The mean average precision as a function of λ 0 is plotted in Figure 7. The mAP values range from 94.16% to 94.29%, yielding a maximum difference of 0.13%, so that it can be concluded that the initial value does not have a large influence on the performance. Note that the best mAP 94.29% is greater than the value reported in Table V. Finally, the effect of the input parameter α on the performance of the similarity measures is investigated. To this aim, the learning algorithm 1 is run with different values of α on the IRIS dataset. The mean average precisions and average losses of algorithms as a function of α varying from 0.1 to 10 are depicted in Figure 8. As can be seen from the graphs, the value α have a larger effect on the performance (up to 0.4%) of the learning algorithm than the initial value λ 0 . As expected, while the loss is decreasing, the corresponding mean average precision is increasing. An interesting point is that local peaks of the loss also correspond to local decreasing of the mean average precision, e.g. log(α) = 1.5, showing the ability of the loss structure to efficiently describe the problem constraints. Here again, one can note that a better performance than the one reported in Table V can be obtained, i.e. greater than 94.18%. Since α is a weight applied on the loss, one can expect that for very small values of α, the loss will be high. On the other hand, when increasing α, the optimal λ value can be reached quickly, but noisy data influences the result.

So far, the number of iterations within the learning algorithm was fixed. The effect of α as a function of the number of iterations is now discussed. The learning algorithm is run for α = 0.1, α = 10 and α = 100. At the end of each iteration, the mean average precision is computed using the actual λ value. Results are given in Figure 9. Naturally, a small α value (α = 0.1) leads to a slow progress rate, since λ i is not very different of λ i-1 . In contrast, when α is large (α = 100), the precision increases faster, but at the price of a worst performance than a medium value (α = 10) later on. Additionally, the more α, the less smooth the curve. A large value of α heavily modifies λ i , resulting in a large difference in terms of performance. In contrast, with a small α value, λ i-1 and λ i slightly differ, which results in a small variation in terms of performance.

V. CONCLUSION

The contribution of this paper lies in the development of generalized fuzzy similarity measures. First, a generic framework of designing similarity measures based on the use of residual implication functions is proposed. This construction presents two main advantages: 1) classical fuzzy similarity measures are retrieved for particular residual functions, 2) verifying if a newly constructed similarity measure satisfies the required properties is facilitated. Then, some new families of parametric similarity measures by using parametric residual implications are presented. An algorithm that allows to learn the parameter of each similarity measure based on relevance degrees is given. Experiments on a number of real datasets show the superiority, which is statistically significant, of the learning algorithm over commonly used similarity functions. The proposed similarity measures can be used in many classification methods, e.g. induction of pattern trees [START_REF] Huang | Pattern trees induction: A new machine learning method[END_REF], hierarchical clustering [START_REF] Mirzaei | A novel hierarchical-clusteringcombination scheme based on fuzzy-similarity relations[END_REF], content-based image retrieval [START_REF] Chan | A region-based fuzzy feature matching approach to content-based image retrieval[END_REF], or ranking image similarities [START_REF] Chechik | Large scale online learning of image similarity through ranking[END_REF].

Among the potential perspectives, a more sophisticated updating scheme, using the passive-aggressive family of algorithms [START_REF] Crammer | Online passive-aggressive algorithms[END_REF] can be developped. In each step, the weight applied on the loss is varying so that the new parameter value better reflects the data. Another interesting perspective is to analyze the fuzzy equivalences between each similarity measures as proposed in [START_REF] Lesot | Order-based equivalence degrees for similarity and distance measures[END_REF] in order to select a particular function within the large variety of measures.

Returning to aggregation operators, it would be interesting to consider ideal samples that characterize the classes, and then adopt a metric that can be learnt with the help of particular aggregation operators [START_REF] Beliakov | Definition of general aggregation operators through similarity relations[END_REF].

Fig. 1 .

 1 Fig. 1. Examples of iso-similarity contours between a reference set A = [ 1 2 ] (denoted by •) and all fuzzy sets F (X), n = 2 where high and low values correspond to white and black colors respectively. The arithmetic mean is used. First row: I ⊤ M , I ⊤ P and I ⊤ L , respectively from left to right. Second row: Hamacher similarity, where λ = 0, 2, 5 from left to right. Third row: Dombi similarity, where λ = 1, 2, 5 from left to right.

Fig. 2 .

 2 Fig. 2. Three fuzzy sets A, B and C where the ranking of similarity measures S(A, B) and S(A, C) is not obvious.

Fig. 3 .

 3 Fig. 3. Three fuzzy sets A, B and C where the similarity of A and B is clearly higher than the similarity of A and C.

4 :

 4 Random selection of a, b and c such that R(a, b) > R(a, c).

Fig. 4 .

 4 Fig. 4. Difference S(A, B)-S(A, C) as a function of the number of update of λ, (see text for details).

Fig. 5 .

 5 Fig. 5. Plot of the 1-dimensional dataset and the membership functions obtained by the fuzzification scheme. Samples of different classes are represented by different colored dots.

Fig. 6 .

 6 Fig. 6. Comparison of each similarity measures for 2 confidence levels α = 0.05 (left) and α = 0.10 (right). In each plot, entries (a, b) are the paired rank differences of similarity measures {i} and {j}, (see text). Black entries correspond to not statistically significant differences under the Nemenyi post-hoc test.

Fig. 7 .

 7 Fig. 7. Mean average precision for the Hamacher based similarity measure S H as a function of the initial value λ 0 .

Fig. 8 .

 8 Fig. 8. Mean average precision (top) and average loss (bottom) for the Hamacher based similarity measure S H as a function of log(α).

Fig. 9 .

 9 Fig. 9. Mean average precision for the Hamacher based similarity measure S H as a function of the number of iterations, for different values of α.

TABLE II ILLUSTRATION

 II OF RANKING CHANGES FOR PARAMETRIC SIMILARITY MEASURES.

	Fuzzy sets	λ	Hamacher similarity
	A, B	2.0	0.505
	A, C	2.0	0.475
	A, B	15.0	0.650
	A, C	15.0	0.762
	Fuzzy sets	λ	Dombi similarity
	A, B	2.0	0.271
	A, C	2.0	0.206
	A, B	0.25	0.686
	A, C	0.25	0.818

TABLE III SIMILARITY

 III MEASURES AND THEIR CORRESPONDING RELATIVE GRADIENT WITH RESPECT TO λ WHEN USING THE ŁUKASIEWICZ IMPLICATION IN THE LOSS FUNCTION.

  . The data contains 3 classes of 50 samples each, where each class is a type or iris plant. Each sample is described by a 4 features: sepal length and width, petal length and width. It is a classical dataset in pattern recognition literature, which is known to have one class linearly separable from the other two, and two classes which are not linearly separable.• the Wine dataset contains 3 classes and 178 samples.Each sample is described by 13 constituents found in each wine. The classes are quite well separated, and samples' classification is not reputed to be a challenging task. • the Pima dataset contains 2 classes of 768 samples, each described by 8 features. The class, 0 or 1, denotes the absence or presence of diabetes pathology. • the Yeast dataset. This is a genetic dataset of phylogenetic profiles for the Yeast genome. It contains 1484 samples, described by 8 attributes. The ten classes of the problem correspond to the various localization sites of the proteins. • the Heart dataset contains 270 samples of 13 attributes each, describing medical data of each patient. The class, 0 or 1, denotes the absence or presence of heart disease. • the Ionosphere dataset is composed of 351 instances. For each sample, there are 17 pulse numbers, each pulse number is described by 2 attributes, giving 34 features. The goal is to discriminate radar returns from the ionosphere.

TABLE IV AVERAGE

 IV PRECISION MEASURE (%) ON IRIS DATASET, FOR THE TOP k SIMILAR SAMPLES (k RANGES FROM 1 TO kMAX = 10).

	k =		1		2	3	4	5		6	7	8	9	10
	S P		92.66	92.66	92	92.33	92.26		92.44	91.71	91.41	90.96	90.53
	S L		94.66	94.66	93.77	92.83	92.53		92.55	92.47	92.41	92.44	92.20
	S M		62.66	62	61.77	62	62.13		60.33	60.19	59.50	58.51	58.20
	S mL		94.66	94.66	94.22	94.16	93.46		92.55	91.61	91.33	90.81	90.60
	S H (λ ⋆ = 58.53) S Y (λ ⋆ = 1.41) S D (λ ⋆ = 0.10) S F (λ ⋆ = 10 7 ) S SS (λ ⋆ = 0.5)	93.33 95.33 95.33 94.66 94.66	94.33 94.66 95 95 94	93.77 94 93.77 94.22 94.44	94.50 93.33 94 93.66 93.66	94.53 92.80 93.46 93.60 93.06		94.66 93 93.44 93.88 93.33	94.47 92.85 93.33 93.61 93.42	94.33 92.83 93.08 93.33 93.33	94.22 92.96 92.96 93.25 93.33	93.66 92.93 92.80 93.26 93.26
								TABLE V
					MEAN AVERAGE PRECISION (%) FOR ALL DATASETS, α = 1.
	S	Iris	Wine	Pima	Yeast	Heart	Ionosphere	Tae	Ecoli	Liver	Newthyroid	Vowel	Breast	Avg. rank
	S P	91.90	93.42	66.45	43.13	75.34	89	48.59	67.80	55.72	91.24	78.27	95.33	7.33
	S L	93.05	95.65	66.39	46.26	75.56	88.97	49.31	74.64	56.05	94.95	82.98	95.87	5.25
	S M	60.73	49.37	52.89	25.21	49.93	70.77	34.26	47.27	49.80	78.82	11.60	80.41	9
	S mL	92.81	92.96	67	45.33	69.06	71.48	48.19	76.45	54.42	91.86	81.96	94.23	7.08
	S H	94.18	94.55	68.39	47.01	75.39	89.64	49.82	76.51	57.73	95.11	85.11	95.97	2.75
	S Y	93.47	95.83	67.32	46.89	75.64	89.66	49.44	74.77	58.37	95.32	83.05	96.02	3.08
	S D	93.72	93.92	68.65	46.79	75.41	92.61	49.54	70.57	57.66	91.97	83.17	95.43	4
	S F	93.85	95.35	68.12	46.77	75.47	89.70	49.43	75.45	55.96	93.65	84.11	95.37	4.08
	S SS	93.65	95.76	67.13	47.12	75.52	89.95	49.46	76.97	59.03	95.44	83.28	95.94	2.41
	OASIS	93.25	94.11	65.41	42.21	73.21	89.58	45.05	73.25	58.05	95.11	82.98	94.37	-


A Banc

B Desbat

D Cavagnat

Ab-Initio calculations of proline vibrations with and without water, consequences on the infrared spectra of proline-rich proteins

Keywords:

Ab-Initio calculations of proline vibrations with and without water, consequences on the infrared spectra of proline-rich proteins

Amélie Banc, Bernard Desbat, D.

Introduction: Proline residues play an important role in the properties and in the structure of globular and fibrous proteins. The absence of amino hydrogen in the chain backbone prevents H-bonded interactions, while the presence of five-membered pyrollidine ring constrains the rotation angle of the C α -N bond. Hence, this residue is often found at bends in the polypeptide chain. Moreover, polyproline can adopt two helical conformations: the form I (PPI) characterized by a right handed helix, in which each peptide bond is in a cis configuration, and the form II (PPII), a left handed helix where the peptide bonds are in the trans configuration 1 . In collagen, prolines seem to pre-organize the individual strands in a PPII conformation to form the typical collagen triple helix [START_REF] Shoulders | [END_REF] . In transmembrane α helices, proline motifs are believed to generate conformational switches at the origin of the translocation of solutes 3 . It has also been suggested that proline residues may play an important role in the folding of membrane proteins 4 . Prolines are also important in prolamins, a class of storage protein of cereal seeds, characterized by their high content of prolines and glutamines (up to 30 mol % and 50 mol % respectively) 5 . These residues are mostly found in repetitive sequences which describe the repetitive domain of these proteins. Spectroscopic studies 6,7,8,9,10,11 indicate that, contrary to the globular non-repetitive domains characterized by classical secondary structures (α-helices, β-sheet), repetitive domains are low folded and subjected to secondary structure changes with the environment (solvent, temperature, pressure…). In infrared studies 12,13 , authors focused on amide I band and assigned the strong band in the 1450 cm -1 region to CH deformation. However, in other studies on bacteriorhodopsin 14,15 and polyproline 16 , authors assigned this strong band to the C-N vibrations of proline. In the 1450 cm -1 region, polyproline infrared spectra exhibit a strong band varying in frequency and intensity according to hydration level 17 and temperature 18 . De Lozé and Josien 17 suggested that these variations were due to formation, by hydrogen bonding, of complex between polyproline and water. Moreover, Caswell and Spiro 16 showed by ultraviolet resonance Raman spectroscopy that polyproline I and polyproline II structures exhibit two different vibration modes (1435 cm -1 and 1465 cm -1 respectively), that they assigned to cis and trans conformations of imide bonds contained in each structure.

In order to establish the 1450 cm -1 band assignment in gliadin infrared spectra, we chose to perform ab initio calculations for a model molecule, C 7 H 13 NO, containing a repeat unit of polyproline (noted I) (Figure 1). This molecule is isomorphous with the structure of cis L proline residue in natural proteins. This simplification is acceptable as we are interested by the vibrations of the proline cycle, which cannot be strongly perturbed by the exchange of the CH 2 group extremities for CH 3 groups. The calculations were firstly done for isolated molecule I. Similar calculations were then done for a complex (I)-H 2 O in order to analyse the hydratation effect observed in previous experimental results obtained for polyproline.

Results:

First of all, we display the infrared spectra of gamma and omega gliadins, two wheat prolamins, in the amide frequency domain (Fig. 2). The two proteins present the characteristic amide I and II absorptions, but also one band around 1450 cm -1 near the CH 2 and CH 3 absorptions. This band is relatively stronger than the bands observed for the other proteins in the same range of frequencies. Moreover, the omega gliadin which contains a higher proportion of prolines displays a lower relative intensity for the amide II and a stronger intensity for the band at 1450 cm -1 . Figure 2. Infrared spectra of lyophilized gamma (continuous line) and omega (dash line) gliadins 11 The evolution of the 1450 cm -1 and of the amide II band peak areas is directly correlated with the proline content in both protein sequences 11 (17.5% of proline in γ-gliadin 19 and 29% in ωgliadin 20 ).

In Table 1, we report the frequencies, infrared intensities and assignment of the B3PW91/cc-pVTZ calculated bands of the isolated molecule I in the amide frequency domain. Figure 3a displays the optimized geometry of this molecule.

Only two strong bands are calculated in this domain. The first one is essentially due to the C=O stretching vibration and is assigned to the amide I band. The second one involves mostly the C-N stretching vibration. In agreement with the structure of the proline residue in the peptidic chain, no amide II band (coupled ν(C-N), δ(N-H) vibrations) is expected, as there is no N-H bond on the tertiary amide group. As a consequence, the decrease of the amide II intensity and the presence of a strong band at lower frequency (ν C-N) in proline rich protein FTIR spectra is directly related with the lower proportion of NH groups in these proteins. This effect can be observed for gliadins but also for collagen proteins, which contain also a high proportion of proline residues.

Frequency

Infrared The B3PW91/cc-pVTZ optimized geometries of molecule (I)-H 2 O complex (figure 3b) shows that the water molecule is hydrogen bonded to the oxygen of the carbonyl. This produces an important modification of the amide frequencies. Indeed, as shown in Table 2, the amide I displays a negative shift of -26 cm -1 and its intensity increases notably in comparison with the vibration of the isolated molecule. Moreover, the ν(C-N) band splits in two close bands with almost the same intensity, but their frequencies are shifted to higher frequencies (+15 and +30 cm -1 ). These results are in good agreement with the experimental spectra obtained for poly-L-proline. With hydration, Wellner et al 12 measured a shift of -26 cm -1 for amide I and +29 cm -1 for the CN vibration mode, and Johnston et al 21 reported shifts of -24 cm -1 and +37 cm -1 respectively. Moreover, when increasing temperature, Swenson et al 18 observed an evolution of frequencies consistent with a decrease of the hydrogen bond strength.

The analysis of the calculated atomic Mulliken charges of molecule (I) reported in table 3 shows that oxygen (3) brings the higher negative charge whereas nitrogen (1) brings the higher positive charge. Thus, the hydrogen bond with water molecule is preferentially done with the oxygen of the carbonyl group. Moreover, the effect of the hydrogen bonded water increases significantly the charges separation between the different atoms and confirms the possibility to get a modification of the amide group vibrational frequencies by the formation of a hydrated system. Table 3. Calculated Mulliken charges for Molecule (I) and complex (I) + H2O

Conclusion:

In summary, the vibration mode observed around 1450 cm -1 for proline rich proteins is characteristic of the coupled ν(C-N), ν(C-C) vibrations of the tertiary amide. The position of this band is sensitive to the hydration state of the oxygen amide group and shifts to higher frequencies when the hydrogen bonds are formed. Moreover, a high proportion of proline residues in one protein produces also a significant decrease of the relative intensity of the amide II band in comparison with the amide I band.

Figure 1 .

 1 Figure 1. Molecule I

  H 20 H 21 (43 %) Table1. B3PW91/cc-pVTZ calculated frequencies (cm -1 ) (* scaled by 0.98) and infrared intensities (KM/Mole) of the isolated molecule I in the 1800-1400 cm -1 domain.

Figure 3 .Table 2 .

 32 Figure 3. B3PW91/cc-pVTZ optimized geometries of isolated molecule (I) (a) and of the molecule (I)-water complex (b)

Proceedings of the National Academy of Sciences of the United States of America 86,9832-9835 (1989).


Olivier Commowick
email: olivier.commowick@irisa.fr

Aymeric Stamm

Romuald Seizeur

Patrick Pérez

Christian Barillot

Sylvain Prima

Nicolas Wiest-Daesslé

Multifiber Deterministic Streamline Tractography of the Corticospinal Tract Based on a New Diffusion Model

In this paper, we build upon a new model, describing the random motion of water molecules in fibrous tissues, to develop a multifiber deterministic tractography algorithm. We apply this algorithm to track the corticospinal tract of the human brain, in both controls and patients with tumors.

Introduction

Tractography of the corticospinal tract (CST) using diffusion-weighted MRI (DW-MRI) is especially challenging, mostly due to the numerous fiber crossings in the corona radiata. When classical diffusion models (e.g. single or multiple tensors) coupled with simple tractography algorithms (e.g. deterministic streamline) are used, these crossings often make it impossible to track the most lateral fibers of the CST [START_REF] Behrens | Probabilistic diffusion tractography with multiple fibre orientations: What can we gain?[END_REF]. These include especially important motor areas such as the hand and the whole face, as shown by the homunculus of Penfield & Rasmussen. The fact that HARDI sequences are prohibitively time-consuming in case of patients with tumors makes it critically important to develop diffusion models and/or tractography algorithms able to track these lateral fibers from clinical (fast) diffusion sequences, having a small number of encoding gradients. We recently proposed a diffusion model that seems to meet these requirements even when using a simple deterministic streamline algorithm [START_REF] Stamm | Diffusion directions imaging[END_REF].

We briefly outline this new model in Section 2.1, the tractography algorithm in Section 2.2, and the pipeline we used to extract the CST in Section 2.3. Finally, we provide tractography results of the left and right CST on the two controls and two patients of the challenge dataset in Section 3.

Methods

Diffusion Modeling

In each voxel, water molecules are assumed to be distributed in several compartments. We first describe how we model the diffusion within a single compartment.

Then, we introduce our multi-compartment model, coined Diffusion Directions Imaging (DDI), and finally we outline how to estimate its parameters.

Single-compartment Model

The diffusion process induces, after a diffusion time τ , a random displacement of water molecules from their initial position x 0 to a random position x = x 0 + √ 2τ w. Assuming a unique direction of diffusion, we propose to model the random variable w as w = u + v, where:

u follows a von Mises & Fisher distribution parametrized by (i) the radius R > 0 of the sphere on which it is defined, (ii) the spherical coordinates (θ, φ) of its mean direction µ and (iii) its concentration parameter κ ≥ 0; v follows a centered Gaussian distribution parametrized by a cylindrically constrained [START_REF] Friman | Uncertainty in white matter fiber tractography[END_REF] covariance matrix D = R 2 κ+1 (I + κµµ ′ ), where I is the identity matrix and {µ, κ, R} are the same parameters that characterize u; u and v are statistically independent.

In essence, (i) µ can be interpreted as the direction of the fibers which constrain the diffusion, (ii) R can be interpreted as the radial displacement along the fiber direction µ and (iii) κ can be interpreted as a measure of anisotropy of the diffusion. The latter can be related to the fractional anisotropy (FA) [START_REF] Basser | Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI[END_REF] as

FA = κ (κ + 1) 2 + 2 -1/2
. We refer the reader to [START_REF] Stamm | Diffusion directions imaging[END_REF] for more details about the motivations of such a parametrization.

The probability density function (pdf) of the molecular displacement xx 0 is then obtained by the convolution of the von Mises & Fisher pdf and the Gaussian pdf, and is parametrized by the four parameters {θ, φ, κ, R}.

Multi-compartment Model Due to its low number of parameters, the singlecompartment model is particularly suited to be encompassed within a multicompartment model, which can account for more than one fiber direction within each voxel. We thereby model the pdf of molecular displacements as a mixture of pdfs having the common parametric form proposed in the previous section.

We assume m compartments associated with m different fiber directions µ i (i = 1, . . . , m). In each compartment, the diffusion is modeled according to the previously described pdf with parameters {µ i , κ i , R i } and mixture weight FA i /m. We also include an additional pdf in the mixture, with weight 1m i=1 FA i /m to account for isotropic diffusion; this pdf follows the general form with κ = 0, so that the unique remaining parameter to estimate is R 0 . Furthermore, we set R 2 i = (κ i + 1)λ, ∀i ∈ 1, m , where λ > 0 is the transverse diffusivity assumed identical in each compartment. For more details about the above described parametrization, we invite the reader to see [START_REF] Stamm | Diffusion directions imaging[END_REF]. Considering m putative fiber directions with this parametrization yields a m-compartment DDI model with 3m + 1 parameters.

Estimation of the DDI Parameters

The theoretical diffusion weighted intensities are the modulus of the Fourier transform of the pdf of molecular displacements which can be analytically derived under the assumption of the DDI model [START_REF] Stamm | Diffusion directions imaging[END_REF]. The 3m + 1 unknown parameters of the m-compartment DDI model are then estimated using a least squares fitting on the raw diffusion weighted intensities, and this optimization is performed using the derivative-free NEWUOA optimization algorithm [START_REF] Powell | The NEWUOA software for unconstrained optimization without derivatives[END_REF]. The different compartments are sorted in decreasing order according to their κ. The model selection is performed according to the procedure described in [START_REF] Stamm | Diffusion directions imaging[END_REF].

Tractography Algorithm

Our goal is to track the fibers linking multiple regions of interest (ROIs). To this end, we developed a deterministic streamline algorithm, which can be viewed as an extension of the original FACT method [START_REF] Mori | Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF], adapted to the DDI model, using a breadth-first-type search.

Starting from one of the ROIs, we define one starting point at each voxel of the ROI. Given one point along its path, we build the main fiber iteratively as follows:

1. If the number of putative fiber directions m = 0, we stop the tracking. 2. If m = 1, we compute FA 1 and the angle α 1 between the input direction and µ 1 . If α 1 < α t and FA 1 > FA t , then we follow the single putative fiber direction µ 1 with a step size of l millimeters. Else, we stop the tracking. 3. If m = 2, we compute FA 1 , FA 2 , the angle α 1 (resp. α 2 ) between the input direction and µ 1 (resp. µ 2 ). If:

α t < α 1 , α 2 : we stop the tracking.

α 1 < α t < α 2 : cf. the case m = 1.

α 2 < α t < α 1 : if FA 2 > FA t , then we follow the direction µ 2 with a step size of l millimeters, else we stop the tracking. -α 1 , α 2 < α t : if FA 2 < FA t , then cf. the case m = 1; else if κ 2 > r × κ 1 then we sort the two fibers in ascending order according to the angles α i . We follow the new direction µ 1 with a step size of l millimeters and we record the second putative fiber direction µ 2 (branch) for future use, as it can be indicative of crossing/kissing/merging/diverging fibers.

Once we have tracked this main fiber, we perform the same tracking from all the possible branching points that we have recorded along its path. Importantly, for these trackings, the stepping rule and stopping criteria are identical as those for the main fiber, but we do not record any possible mixed fiber configuration along these secondary paths, for which we only follow the main direction µ 1 at each step. We then lead the same tracking from the other ROIs. We only keep the tracts going through all the seeding ROIs for further analysis. In practice, we choose the parameters l = 1, r = 0.8, α t = 60 degrees (maximal angle between two successive directions along the fiber) and FA t = 0.5 (minimal FA along the fiber). During the tracking, when a point is not on the grid of the DW-MR images, we compute the DDI model using a trilinear interpolation directly on the model parameters.

Tractography Pipeline Applied to the Challenge Datasets

We utilize the following processing pipeline to extract the CST for the challenge datasets. An expert neuroanatomist (Romuald Seizeur) delineated two ROIs on each side (left and right) of the original T1-weighted images. One is located in the posterior limb of the internal capsule and the other in the superior part of the mesencephalon. In addition, since crossing fiber tracts such as the association or commissural fiber tracts may be considered as bifurcations of the CST by the tractography algorithm, the same expert also delineated regions through which the tracts are not allowed. On the patients' datasets, the same ROIs were delineated with some modifications to account for the deformations caused by the tumor (e.g. in patient 1 the internal capsule on the left side is compressed by the infiltrating tumor).

All datasets were then processed in three steps:

1. Diffusion-weighted MRI denoising: DW-MRI is subject to random noise yielding measures that are different from their real values, and thus biasing the subsequently estimated diffusion models. We filtered the diffusionweighted MR images with the Rician-adapted Non-Local Means filter [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio MRI: applications to DT-MRI[END_REF], which has been shown to efficiently denoise such images while preserving fine anatomical structures. In particular, this filter has also been shown to preserve the angular resolution of q-ball ODF models estimated from HARDI data [START_REF] Descoteaux | Impact of Rician Adapted Non-Local Means Filtering on HARDI[END_REF]. 2. ROI alignment on B0 images: we registered the ROIs on the B0 images according to the following steps:

global affine registration of the T1-weighted images to the B0 images [START_REF] Ourselin | Block matching: A general framework to improve robustness of rigid registration of medical images[END_REF]; -cropping of the affine-registered T1-weighted images using the mask of the B0 images; -constrained non-rigid registration [START_REF] Garcia | A Robust and Efficient Block Matching Framework for Non Linear Registration of Thoracic CT Images[END_REF] of the masked T1-weighted images to the B0 images; -application of the obtained transformations to the ROIs. 3. Extraction of the left and right CST using the aligned ROIs and the DDI estimated from the DWI.

Results

The MICCAI DTI tractography challenge consists of two groups of data: two healthy subjects acquired using a multiple b-values scheme and two patients acquired using a more standard acquisition protocol. For each of the two healthy subjects, DW-MRI data were acquired repeatedly (ten repetitions) so that the reproducibility of the tractography method may be evaluated. On the other hand, the tumors of the two patients were delineated so that we could produce combined views to help the surgeon, for example when planning a tumor removal surgery.

Tractography on Healthy Subjects

For each of the repetitions of each subject, the left and right CST were computed utilizing the aligned ROIs (see Section 2.3). We display the tractography of one volume of each of the two subjects in Fig. 1. This figure illustrates that we are able to cover the full extent of the CST for these healthy subjects, from the face area to the hand area to the medial part of the CST. This demonstrates that our diffusion model enables the tractography algorithm to follow bifurcations in the white matter fiber tracts. For each subject, estimating the DDI models from the DWI took approximately 40 minutes (singlethreaded), while the tractography of the left and right CST took less than ten minutes.

Tractography on Patients

We also report results for the CST of each patient on the side of the tumor (all fiber bundles are provided in supplemental material). The images processed here had a larger resolution than the healthy subjects and the estimation of the DDI models from the DWI took approximately 3 hours (single-threaded). However, this step of the tractography pipeline may be computed offline leaving only the tractography to perform online (about five minutes for each CST) when performing the tractography for surgery planning.

We focus on the qualitative evaluation of the obtained tracks and their closeness to the tumor (it should be noted that the tumor regions were not used in any way to constrain the tractography algorithm), and on providing the neurosurgeon with helpful views for neurosurgery planning. To this end, we present in Fig. 2 views (made using the MedINRIA software [START_REF] Toussaint | MedINRIA: Medical Image Navigation and Research Tool by INRIA[END_REF]) combining the fiber tracts, the tumor ROIs provided by the organizers, all on top of the patient's T1 image. More illustrations are available in supplemental material.

We can observe on this figure that, although the tumor delineation was not used in the algorithm, no fibers are going through the tumor area on patient 1. Instead, the CST is going through the most central part of the edema (see images (b,c) on Fig. 2). This indicates that the fibers were pushed by the tumor mass effect, which is a valuable indication when planning the surgery. Overall, this figure demonstrates the close proximity of the tumors and of the CST for both patients. This is an important insight as the neurosurgeon will be able to plan the tumor removal in the optimal way, so as to minimize the possible handicap for the patient after surgery.

Fig. 1 .

 1 Fig. 1. Illustration of CST on Healthy Subjects. Representative examples of obtained fiber tracts for healthy subject 1 (a) and 2 (b) (T1 images are in radiological conventions, i.e. the left hemisphere is on the right side of the image). To see the full extent of fiber spreading, all 3D fiber tracts are displayed. This explains why they may not seem to match exactly the background T1 image.

Fig. 2 .

 2 Fig. 2. Illustration of CST on Patients. Combined views of the obtained fiber tracts for patients 1 (a,b,c) and 2 (d,e,f), illustrating the proximity of the tumor to crucial motor pathways. Surfaces for the first patient correspond to the necrotic part of the tumor (red), the active part of the tumor (yellow) and the edema (orange). For the second patient, each ROI corresponds to a specific tumor. Images (a,d) show overall 3D views and (b,c,e,f) show the tracts and regions of interest going through a specific 2D slice to better illustrate their proximity.

Acknowledgments This investigation was partly supported by the ARSEP (French MS society, Fondation pour l'Aide à la Recherche sur la Sclérose en Plaques).