Keywords: Diffusion MRI, White Matter Fiber Tracts, Gaussian Process, Tract Probability Maps, Clustering IRM de Diffusion, Imagerie par Tenseur de

Regroupement et Carte de Probabilités des Faisceaux de Fibres de la Matière Blanche : Application du concept de Processus Gaussien à la description des Fibres

Résumé : Au vu de l'intérêt clinique croissant pour la tractographie des fibres en imagerie du tenseur de diffusion, la demande d'un cadre mathématique objectif pour l'analyse quantitative des fibres de la matière blanche, incorporant leur signification physique sousjacente, s'est faite de plus en plus pressante. Cet article présente un tel cadre, facilitant les opérations mathématiques entre les faisceaux grace à un produit interne (basé sur la notion de processus gaussien) entre les fibres qui définissent alors un espace métrique. Cette métrique facilite la combinaison de faisceaux de fibres, rendant les calculs d'appartenance à un faisceau ou de similarités entre faisceaux simples. À partir de ces notions, nous avons défini une méthode automatique et non-supervisée de regroupement de faisceaux qui se base sur un atlas, et qui ne nécessite ni initialisation manuelle, ni connaissance a priori du nombre de faisceaux. On peut dès lors faire une analyse quantitative inter-sujets sur le volume des faisceaux regroupés, en évitant ainsi de devoir paramétriser les fibres, de recourrir à l'enveloppe du faisceau, ou à la fibre moyenne, contrairement aux études précédentes. Les opérations mathématiques sont montrées sur des données synthétiques, puis suit une étude de l'applicabilité de cette méthode de regroupement non supervisé, sur un jeu de 21 sujets.

Introduction

Diffusion MRI non invasively recovers the in vivo effective diffusion of water molecules in biological tissues. This information characterizes tissue micro-structure and its architectural organization [START_REF] Basser | Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI[END_REF] by modeling the local anisotropy of the diffusion process of water molecules, providing unique biologically and clinically relevant information not available from other imaging modalities. Once the diffusion information has been recovered within each voxel, it can be synthesized in the form of a diffusion tensor [START_REF] Basser | MR diffusion tensor spectroscopy and imaging[END_REF]. Brain connectivity can then be assessed by assembling the tensors into tracts using tractography methods [START_REF] Mori | Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF][START_REF] Koch | An Investigation of Functional and Anatomical Connectivity Using Magnetic Resonance Imaging[END_REF][START_REF] Descoteaux | Deterministic and probabilistic tractography based on complex fiber orientation distributions[END_REF]. Among these methods, streamline tractography [START_REF] Mori | Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF] recovers white matter fiber tracts from a seed voxel by following the principal direction of the diffusion tensor. Using this technique, white matter fiber tracts are represented as points sampled from a three-dimensional curve. Finally, these fibers can then be grouped into fiber bundles based on anatomic knowledge (O'Donnell and Westin, 2007;Maddah et al., 2008a). The Cortico Spinal Tract (CST) or the Corpus Callosum (CC) are prominent examples of the latter.

In this paper, we address the important problem of developing a mathematical framework for the quantitative analysis of fiber bundles, which has become a very active research area, with the aim of facilitating subsequent clustering and group-based statistical analysis on the bundles. The effects of pathology on tracts is evident, as in the case of brain tumors that grossly displace tracts, or can be subtle, as in the case of neuropsychiatric disorders, such as schizophrenia, which can manifest as more diffuse abnormalities in connectivity [START_REF] Kubicki | A review of diffusion tensor imaging studies in schizophrenia[END_REF][START_REF] Ciccarelli | Diffusionbased tractography in neurological disorders: concepts, applications, and future developments[END_REF]. We propose a framework for tract analysis based on a novel metric between bundles that serves as a probabilistic measure of inclusion of a fiber tract into a bundle. This framework is then used to develop a method for automated clustering of bundles, which are now quantifiable on the basis of membership and ready for statistics based on the distances between bundles. Once obtained, the clusters can be mapped to tract probability maps [START_REF] Hua | Tract probability maps in stereotaxic spaces: Analyses of white matter anatomy and tract-specific quantification[END_REF] enabling tract-based statistics on the cerebral white matter.

Materials and Methods

Prior work

The clustering of different fiber tracts into an anatomically coherent bundle, like the CC or the CST, is a challenging task for several reasons. In the first place, as seen in fig. 1, axons composing a bundle can diverge from it connecting cortical and subcortical areas. This renders approaches that quantify similarity among white matter fibers using the whole fibers instead of analyzing partial overlaps like shape statistics or rigid transformations [START_REF] Veltkamp | Shape matching: Similarity measures and algorithms[END_REF] unsuited for the clustering task. Take for instance the cingulum bundle, whose constituent fibers only partially overlap among themselves, with many diverging to innervate the cortex, as shown in fig. 1. These divergent fibers can have quite different shapes, calling

INRIA S A
Anatomical region of the cingulum.

Axons actually enter and leave.

Figure 2-6: Axons enter and leave some fiber tracts along their course, as seen in this example of the cingulum bundle. Images courtesy of James Fallon, Professor of Anatomy and Neurobiology at the University of California, Irvine.

have been reported in schizophrenia [63], and asymmetry of the uncinate fasciculus across hemispheres was found in a postmortem dissection study [48].

The Inferior Occipitofrontal Fasciculus

The inferior occipitofrontal fasciculus (Figure 2-5) connects occipital and frontal lobes.

It runs from the frontal lobe, through the temporal lobe (passing below the insula), to the occipital lobe [84].

The Inferior Longitudinal Fasciculus

The into question the utility of shape-based metrics for subsequent tensor statistics. Moreover, even if current streamline tractography techniques are reproducible [START_REF] Wakana | Reproducibility of quantitative tractography methods applied to cerebral white matter[END_REF], fiber tracts obtained through tractography do not recover the whole underlying axonal trajectory [START_REF] Lenglet | Mathematical methods for diffusion mri processing[END_REF]. This behaviour can be clearly observed in complex bundle configurations like crossings or fannings [START_REF] Savadjiev | Labeling of ambiguous subvoxel fibre bundle configurations in high angular resolution diffusion MRI[END_REF]. An example of this is the crossing area of the CC and the CST. Thus, despite innovative techniques that are more sensitive to fiber tracts being used in order to model the underlying water diffusion in a voxel and produce better local approximations to axonal distribution [START_REF] Qazi | Resolving crossings in the corticospinal tract by two-tensor streamline tractography: method and clinical assessment using fMRI[END_REF][START_REF] Descoteaux | Deterministic and probabilistic tractography based on complex fiber orientation distributions[END_REF][START_REF] Wiegell | Fiber Crossing in Human Brain Depicted with Diffusion Tensor MR Imaging 1[END_REF], results are far from being reliable after the streamline tracking procedure has traversed a region with complex bundle configuration.

In fiber bundle identification (clustering tracts into a bundle), a common workaround for this problem consists of performing a dense whole brain tractography by densely seeding all over the brain. This procedure produces fiber tract sections that later could be grouped by clustering techniques (O'Donnell and Westin, 2007). This calls for the need for a similarity metric, that can quantify the closeness of two fibers or the degree to which a tract belongs to a bundle, which are both challenging problems and critical in automatic bundle integration. Quantifying fiber similarity is a fundamental part of fiber clustering that has been addressed in different ways. Recent works [START_REF] Batchelor | Quantification of the shape of fiber tracts[END_REF][START_REF] Corouge | Fiber tract-oriented statistics for quantitative diffusion tensor mri analysis[END_REF][START_REF] Leemans | Multiscale white matter fiber tract coregistration: a new feature-based approach to align diffusion tensor data[END_REF] quantify fiber similarity with different flavors of shape statistics. However, partial overlapping of fibers is not taken into account as a similarity feature. Thus, the previous approaches are unsuited for automatic classification of fibers in the brain. There is a separate set of works [START_REF] Ding | Classification and quantification of neuronal fiber pathways using diffusion tensor MRI[END_REF]O'Donnell and Westin, 2007;[START_REF] Wassermann | Simultaneous manifold learning and clustering: Grouping white matter fiber tracts using a volumetric white matter atlas[END_REF]Maddah et al., 2008a), which uses different clustering algorithms based on the Hausdorff or Chamfer distances among the sequence of points parametrizing each fiber tract. This family of similarity metrics deals with sets of points instead of curves, hence they discard continuity or directionality information. Moreover, similarity tends to decrease very fast in cases of partial overlapping, failing to include fibers diverging from the bundle in the correspondent cluster. In particular, [START_REF] Ding | Classification and quantification of neuronal fiber pathways using diffusion tensor MRI[END_REF] only analyzes fibers whose seed points are spatially close together. This is not suited for a whole brain analysis. In this case different fiber seed points from the same bundle might have been scattered all over the white matter. The latter technique is frequently used in order to overcome limitations of streamline tracking in regions with complex bundle configurations (O'Donnell and Westin, 2007). Manifold learning techniques are used by O'Donnell and Westin (2007) and [START_REF] Wassermann | Simultaneous manifold learning and clustering: Grouping white matter fiber tracts using a volumetric white matter atlas[END_REF] to generalize these type of distances from small sets of similar fibers to a bigger more diverse set of fibers. This approach embeds the fibers into Euclidean or topological spaces that can be handled more easily. This strategy has proved to be useful for single individuals but lacks the necessary parameter stability needed for group studies. O' Donnell and Westin (2007) start by generating an atlas, fibers from new subjects are then classified according to this atlas. Even though these automatically grouped bundles are anatomically coherent, the process to generate the atlas requires heavy user interaction and fine parameter tuning. The level of manual interaction needed renders the approach difficult to reproduce. Wassermann and Deriche (2008) use a publicly available anatomical atlas in conjunction with the fiber similarity metric. This work requires a smaller number of parameters, nevertheless situations of partial fiber overlapping generate non-anatomically coherent bundles. Maddah et al. (2008a) enhances the Hausdorff similarity with Mahalanobis distance between fiber points. In order to handle partial overlapping, an ad-hoc penalty term is added to this distance. This approach requires user initialization, by selecting a fiber which is known to be in the desired bundle. Their subsequent work (Maddah et al., 2008b) incorporated atlas information to increase accuracy, however an initial fiber representing each bundle is still required. From all the presented approaches, only O' Donnell and Westin (2007) succeed in the task of semi-automated classification of the whole ensemble of white matter fibers, however this is achieved with a great deal of user interaction and parameters tuning.

Once bundles have been found, quantitative analysis can be performed. This analysis is useful to monitor pathological conditions [START_REF] Ciccarelli | Diffusionbased tractography in neurological disorders: concepts, applications, and future developments[END_REF][START_REF] Kubicki | A review of diffusion tensor imaging studies in schizophrenia[END_REF]. Most of the works performing bundle statistics [START_REF] Goodlett | Group statistics of dti fiber bundles using spatial functions of tensor measures[END_REF][START_REF] Hua | Tract probability maps in stereotaxic spaces: Analyses of white matter anatomy and tract-specific quantification[END_REF]O'Donnell et al., 2007;Maddah et al., 2008a) rely on the use of medial representations for bundles. These representations are only appropriate for bundles which can be modeled as convex envelopes. Thus, their methodology is not entirely appropriate at the extremes of the bundles, an area where the axons fan-out inervating cortical or subcortical structures. Other works [START_REF] Oh | Tractography-guided statistics (tgis) in diffusion tensor imaging for the detection of gender difference of fiber integrity in the midsagittal and parasagittal corpora callosa[END_REF] use a mesh approach over ROIs and perform statistics on the surface. However, automated transition from a set of fibers to the mesh is unclear. There is recent evidence that tract probability maps can be used in order to perform bundle-oriented statistics in diffusion MRI [START_REF] Hua | Tract probability maps in stereotaxic spaces: Analyses of white matter anatomy and tract-specific quantification[END_REF] and histological images [START_REF] Bürgel | White matter fiber tracts of the human brain: Three-dimensional mapping at microscopic resolution, topography and intersubject variability[END_REF]. However, in these two approaches the process to obtain the tract probability maps must be performed manually by experts. Overall, statistical models of white matter bundles which do not rely on medial representations are scarce. Moreover, for models which are more appropriate for a wider spectra of bundles, the transition from automatically obtained bundles to these statistical models is not straightforward.

INRIA

It is the goal of this work to introduce a novel mathematical framework to model the combination of fiber tracts and bundles. Our model includes diffusion information and relates the bundles with a ROI in the volume, mapping every voxel to its probability of belonging to the bundles, the tract probability map. Our framework provides an inner product space for fiber bundles and fiber tracts, which are considered as single-fiber bundles. Inherently, it provides a metric between these elements. Moreover, linear combinations of fiber bundles seamlessly generate new bundle configurations and allow for the volume-based statistics of fiber bundles. In addition, the framework handles cases of partial fiber overlap, in fiber combinations, naturally. Furthermore, the previously mentioned characteristics of our framework facilitate statistical analysis and classification/clustering tasks. Finally, we present a clustering application based on our mathematical framework and on anatomical information in the shape of a volumetric atlas. The output of this application is a set of automatically obtained white matter bundles like the Arcuate Fasciculus or the Cingulum. For each bundle, we are also able produce a ROI which maps every voxel to its probability of belonging to the bundle, referred in previous work as tract probability map [START_REF] Hua | Tract probability maps in stereotaxic spaces: Analyses of white matter anatomy and tract-specific quantification[END_REF][START_REF] Bürgel | White matter fiber tracts of the human brain: Three-dimensional mapping at microscopic resolution, topography and intersubject variability[END_REF]. This map is an appropriate tool to perform bundle-based statistics on the cerebral white matter. We validate this clustering algorithm by applying it to 21 healthy subjects and then performing statistical analysis on the results using our mathematical framework.

A Mathematical framework for white matter fibers and bundles

In this section, we develop our mathematical framework for the quantitative analysis of fiber bundles with the aim of facilitating subsequent clustering and group-based statistical analysis on the bundles. Within this framework, examples of operations such as the combination of fibers into bundles, similarity measurements between fibers, as well as tract probability maps computation will be easily provided. To this end, we first need to model each fiber as a stochastic process. This is done using the well known Gaussian Process framework [START_REF] Mackay | Introduction to gaussian processes[END_REF] which has now a long history in the field of statistics. Indeed, a Gaussian process can be seen as a generalization of the classical Gaussian probability distribution to describe properties of functions and not only properties of random variables such as scalars or vectors. A stochastic process will be defined to model each fiber as a function. This will allow us to govern the properties of our fibers/functions while subsequently simplifying a certain number of sophisticated operations between the fibers. Examples of linear combination of Gaussian Processes of several fibers resulting in the Gaussian Process of a bundle are shown in fig. 3. It is important to point out that such framework provides also an adequate computational tractability. It allows us to measure bundle similarity with partial fiber overlap without relying on point correspondences as shown in fig. 4. Furthermore, the linear combination and inner product operations span an inner product space, which is an appropriate framework to perform statistical analyses such as clustering on fiber bundles. Finally, as we model bundles with stochastic processes, producing a tract probability map for each bundle is a natural operation. Thus, the mathematical framework we develop in this section provides all the right and necessary tools to perform a quantitative statistical analysis of white matter fiber bundles. We model each fiber tract F as a blurred indicator function y F : p∈R 3 →R. This blurred indicator function has a maximal level set which corresponds to F. Moreover, y F is blurred in accordance with the information provided by the underlying diffusion tensor field, i.e. along the fiber direction and not across it. This is illustrated on fig. 2. As we show in appendix A, Gaussian Processes (GP) provide the right framework to integrate spatial and diffusion information for y F . More precisely, the blurred indicator function y F (p) is modeled by the GP

y F (p)∼GP(y * F (p), c F (p, p )), (1)
where the mean function y * F (p) and covariance function c F (p, p ) are the parameters of this stochastic process. These two functions are inferred from the tractography of each fiber. That is, from the sequence of points f = {f 1 , . . . , f |f | } estimated by tractography of the anatomical bundle F and from the corresponding sampling on its tensor field Σ(f 1 ), . . . , Σ(f |f | ). The inference process of the GP corresponding to a white matter fiber is developed in appendix A.1. Using our framework, we are able to combine N fibers into a bundle with a simple average operation in our inner product space. The GP which corresponds to the indicator function of a fiber bundle B, is obtained through the mean Gaussian Processes of single fibers or smaller bundles composing it,

y B (p) = 1 N N i=1 y Fi (p)∼GP (y * B (p); c B (p, p )) , (2)
where the value of our inner product operation F, F , defined in eq. 3, and of our inner product normalized by its natural norm, F = F, F . Inner product quantifies the overlapping of blurred indicator functions. A larger inner product means that fibers are more similar and relates to the volume of the overlapping. The inner product normalized by its norm quantifies similarity ranging from 0 when overlapping is null to 1 when the two fibers are identical. The compared fibers have been extracted from different anatomical tracts of a diffusion MRI image, the Frontal Forceps (FF), the Uncinate Fasciculus (UNC) and the Cingulate Cortex section of the Cingulum (CgC).

y * B (p) = 1 N N i=1 y * Fi (p) and c B (p, p ) = 1 N 2 N i=1 c Fi (p, p ). INRIA (a) F 1 (b) F 2 (c) F 3 (d) F 4 (e) P i F i 4
An example of this is shown in fig. 3, where indicator functions for four fiber tracts and the obtained function for the bundle combining them can be seen. Within the proposed framework, we measure similarity between bundles by quantifying their overlapping. This is done with an inner product operation developed in detail in appendix A.3. If the mean function for a bundle F is square integrable and has finite support, the inner product between two bundles F and F is defined as

F, F := R 3 y * F (p)y * F (p)dp, (3)
along with its induced norm

F 2 := F, F . (4)
Moreover, we can define a similarity measure, bounded by 1 when F and F are exactly the same and by 0 when there is no intersection, as the inner product normalized by its induced norm,

F, F N := F, F F F . ( 5
)
Examples of these two similarity measures are shown in fig. 4. We mentioned previously the importance for bundle statistics of the tract probability map, the probability that a point p in R 3 is contained in a bundle B, P {p ∈ F}. As we show in appendix A.4, having calculated the GP for a bundle B, y B (•), by means of eq. 2, we can define the tract probability map from the parameters of the GP as

P {p∈F} ∝ 1 2 π(h 2 + σ 2 B (p) , (6)
where h is a parameter that diffuses the tract probability map in space and σ 2 B (p) is calculated from the parameters of y B (•). In order to illustrate the tract probability map calculated from a GP for a bundle, color-coded surfaces and a probability map over an FA image for anatomical bundles are shown in fig. 5.

Up to this point, we have introduced our GP-based framework for white matter fiber bundles and its three main operations: combination of fibers into a bundle, similarity quantification and calculation of the tract probability map. In the following sections we use these tools in order to perform automatic white matter bundle identification by means of white matter fiber clustering and we assess the quality of this clustering by performing inter-subject statistical analysis of white matter fiber bundles using our framework.

Clustering Algorithm

Taking advantage of the mathematical framework for fiber bundles presented in section 2.2, we propose a stochastic process-based agglomerative clustering algorithm. This algorithm is executed over a set of densely sampled fibers from the whole white matter, a full brain tractography. Once executed, our algorithm generates a dendrogram, a tree structure where each joint is a candidate cluster, this is illustrated in fig. 6. Then, this dendrogram can be interactively explored in order to choose the desired granularity of the clustering without reprocessing the data. Furthermore, being a hierarchical agglomerative algorithm it has several desirable properties. To begin with, convergence is guaranteed by the finite number of elements to cluster. Next, the number of clusters does not have to be known a priori. Furthermore, outliers are incorporated to clusters during the late stages of the process if at all, thus being easily distinguishable from clustered bundles [START_REF] Jain | Data clustering: a review[END_REF]. All of these characteristics make our clustering algorithm effective and robust in order to classify white matter fibers from a full brain tractography into anatomically coherent bundles.

Our clustering algorithm applied to a full brain tractography is as follows,

Step 1. Given a full brain tractography F = {F i }, with 1 ≤ i ≤ |F |, calculate the set of stochastic processes representing each fiber Y = {y Fi (p)}.

Step 2. Initialize the clustering as the set of single fiber bundles B = {B i }, where B i is a set (bundle) formed by a single fiber, F i , and

1 ≤ i ≤ |F |.
Step 3. Initialize the set of edges in the dendrogram: T = ∅

Step 4. While there is a pair of different bundles B,B in B, s.t. B, B > 0

Step 4.1. Select two different bundles B, B such that B, B = max C,C ∈B C, C .

Step 4.2. Remove the bundles B and B from B and add the bundle {B + B }.

Step 4.3. Add the edge (B, B ) to the dendrogram T .

The output of this algorithm is a dendrogram T , more precisely set of trees where each edge represents the joining of two bundles such as the one shown in fig. 6.

A main advantage of our framework within this clustering algorithm is that the most important operation in step 4, the inner product among bundles described in section 2.2, ... is fast and simple to compute as we show in appendix A.4. Thus, once we have calculated the matrix of inner products for every pair of fibers in the full brain tractography F, the algorithm works by simply performing linear operations on the rows of this matrix.

Tract Querying: Automatic cluster selection based on anatomical knowledge

Once our clustering algorithm generates the dendrogram for a full brain tractography, section 2.3, the main problem is how to select the joints in the dendrogram so that they are anatomically correct clusters. In order to do this, we introduce a query system based on volumetric information, the result of this query will be a cluster selected from joints of the dendrogram. In this work we use a publicly available atlas which has a parcellation of the brain gyri on the grey and white matter [START_REF] Wakana | Fiber Tract-based Atlas of Human White Matter Anatomy[END_REF] as anatomically-aware volumetric information. Then, an anatomical query is defined by a set of grey or white matter regions that the tract must traverse, for example, the Inferio Fronto Occipital tract must connect the inferio-frontal gyrus and the medial-occipital gyrus. Images of the white matter atlas parcellation and queries for various tracts are shown in fig. 7. After setting an anatomical query Q traversing K labeled regions on the atlas, Q = {r 1 , . . . , r K }, we use the tract probability map of each bundle, eq. 6, to select the bundle B on the dendrogram T with 3 Experiments

INRIA

Subjects, Imaging and Data Processing

Whole-brain DWI datasets were acquired from 21 healthy volunteers (30.05 +-7.05 years, 10 Male ) on a Siemens Trio 3T scanner with 1.71x1.71mm 2 in-plane resolution, 2mm thick slices, six unweighted images and 64 diffusion weighted images (b=1000s/mm 2 ) acquired with non-collinear diffusion sensitizing gradients. DTI images for each subject were computed and deformably registered, using DTI-DROID [START_REF] Yang | Diffusion Tensor Image Registration Using Tensor Geometry and Orientation Features[END_REF], to a DTI atlas [START_REF] Wakana | Fiber Tract-based Atlas of Human White Matter Anatomy[END_REF]. Full brain tractography was performed following (O'Donnell and Westin, 2007), streamline tractography was performed by seeding in sub-voxel resolution by taking every voxel with linear anisotropy higher Using our new mathematical framework for white matter bundles, we develop a clustering scheme to perform automatic extraction of white matter fiber bundles across subjects. Firstly, using our similarity measure we build an algorithm which receives a full brain tractography as an input and outputs a dendrogram, Figure 5. The dendrogram, a precomputed tree structure where each joint is a candidate cluster, is a versatile tool to represent clustering results. Among other advantages, it allows interactive exploration of the clustering results at different scales and it is resilient to outliers. Moreover, clusters can be selected at different scales using different postprocessing methodologies. In particular, we take advantage of this feature using the tract prob- DTI images for each subject were computed and deformably registered, using DTI-DROID [START_REF] Yang | Diffusion Tensor Image Registration Using Tensor Geometry and Orientation Features[END_REF], to a DTI atlas [START_REF] Wakana | Fiber Tract-based Atlas of Human White Matter Anatomy[END_REF]. Full brain tractography was performed following (O'Donnell and Westin, 2007), streamline tractography was performed by seeding in sub-voxel resolution by taking every voxel with linear anisotropy higher than .3, dividing it into .25x.25x.25mm 3 sub-voxels and seeding from each sub-voxel. In average, 10.000 fiber tracts where obtained for each subject. Then, the previously presented clustering algorithm based on our mathematical framework was applied to every subject individually. In order to extract major white matter tracts on every subject individually we performed set of queries, shown in Table 1, over the dendrogram obtained from the clustering of each subject using the tract querying algorithm we introduced.

Results

We show the results of our clustering and tract querying algorithms for 4 different subjects selected from our 21 subject database in Figures 7 to 10. Then, we evaluated cluster coherence across subjects. Firstly, calculation of the GP representation of the population-averaged bundle across the 21 subjects for each automatically extracted tract: For each queried white matter tract, the GP corresponding to the populationaveraged bundle across subjects was calculated using Equation 2. We show the tract probability map for each population-averaged bundle, calculated with Equation 6, in Figures 11 and12. These maps are in atlas space, over fractional anisotropy images. Finally, in order to provide quantitative evaluation, for each bundle extracted from each subject we quantified its similarity with respect to the corresponding populationaveraged bundle using the normalized inner product of our mathematical framework, Equation 2.2. The result is plotted in Figure 13. In this plot, the boxes span between the second and third quartiles of the similarity with the population-averaged bundle, the red bar is the median similarity. Moreover, the whiskers indicate the bundles whose 10 The procedure used to cluster white matter fibers into anatomical bundles and produce the tract probability maps for each bundle: We registered 21 Diffusion Tensor MR images using [START_REF] Yang | Diffusion Tensor Image Registration Using Tensor Geometry and Orientation Features[END_REF]. We performed full brain tractography obtaining around 10.000 fibers per brain. We produced the Gaussian Process representation for each fiber as we describe in section 2.2. We identified anatomical bundles like the arcuate or the uncinate fasciculus by applying our clustering and tract querying algorithms to each subject individually as we describe in section 2.3. Finally we produced a population-averaged Gaussian Process for each identified bundle and the corresponding tract probability map with the methodology described in section 2.2 than .3, dividing it into .25x.25x.25mm 3 sub-voxels and seeding from each sub-voxel. In average, 10.000 fiber tracts where obtained for each subject. Then, the previously presented clustering algorithm based on our mathematical framework was applied to every subject individually. In order to extract major white matter tracts on every subject individually we performed set of queries, shown in table 1, over the dendrogram obtained from the clustering of each subject using the tract querying algorithm we introduced.

( a ) C o r t i c o S p i n a l T r a c t L e f t ( b ) C o r t i c o S p i n a l T r a c t R i g h t ( c ) U n c i n a t e F a s c i c u l u s L e f t ( d ) U n c i n a t e F a s c i c u l u s R i g h t ( e ) S u p e r i o r L o n g i t u d i n a l F a s c i c u l u s L e f t ( f ) S u p e r i o r L o n g i t u d i n a l F a s c i c u l u s R i g h t ( g ) C i n g u l u m L e f t ( h ) C i n g u l u m R i g h t ( i ) I n f e r i o F r o n t o O c c i p i t a l L e f t ( j ) I n f e r i o F r o n t o O c c i p i t a l R i g h t F i g u

Results

We show the results of our clustering and tract querying algorithms for 4 different subjects selected from our 21 subject database in fig.s 9 to 12. Then, we evaluated cluster coherence across subjects. Firstly, calculation of the GP representation of the population-averaged bundle across the 21 subjects for each automatically extracted tract: For each queried white matter tract, the GP corresponding to the population-averaged bundle across subjects was calculated using eq. 2. We show the tract probability map for each population-averaged bundle, calculated with eq. 6, in fig.s 13 and 14. These maps are in atlas space, over fractional anisotropy images. Finally, in order to provide quantitative evaluation, for each bundle extracted from each subject we quantified its similarity with respect to the corresponding population-averaged bundle using the normalized inner product of our mathematical framework, eq. 2.2. The result is plotted in fig. 15. In this plot, the boxes span between the second and third quartiles of the similarity with the population-averaged bundle, the red bar is the median similarity. Moreover, the whiskers indicate the bundles whose similarity value with the population-averaged bundle is the smallest and largest within within 1.5 times the interquartile distance of the population-averaged similarity and the '+' symbols indicate outliers. 1. The subjects where previously registered to a white matter fiber atlas and full brain tractography was performed by densely seeding on the whole cerebral white matter. Four subjects are shown.

Discussion

In recent years tractography has become a popular means to perform white matter studies through diffusion MRI [START_REF] Ciccarelli | Diffusionbased tractography in neurological disorders: concepts, applications, and future developments[END_REF]. Tractography results are visually appealing and recent developments produced means to perform analysis of diffusion-derived measures Maddah et al., 2008a) and tract probability maps have been previously reported [START_REF] Hua | Tract probability maps in stereotaxic spaces: Analyses of white matter anatomy and tract-specific quantification[END_REF]. However, clustering approaches were not suited for large-scale clinical studies or required a great deal of parameter tuning. With regard to tract probability maps, they were calculated by averaging binary masks obtained by manually extracted white matter fibers from tractography or by chemical staining. This paper proposes a mathematical framework that provides the necessary tools to perform automatic clustering and identification of white matter fibers and subsequent statistics on them. Thus, automatic clustering of white matter fibers into bundles which is consistent among individuals is a fundamental tool for clinical studies. In this work we provide the tools to solve these issues. Firstly, we develop a mathematical framework for handling white matter bundles which includes spatial and diffusion information and provides the grounds for their statistical analysis. Secondly, we use this framework to develop a clustering algorithm and automatically obtain white matter structures among subjects, which is then applied to 21 subjects.

Our mathematical framework for white matter bundles, sets the grounds for statistical analysis of white matter fiber bundles including applications like clustering and tract-based quantification of scalar quantities among others. It provides three important operations for white matter bundle statistics: 1) an operation to merge white matter fibers into a bundle by means of a linear combination within our framework, see Fig. fig. 3; 2) an operation to quantify bundle similarity based on their overlap in space (see Fig. fig. 4, quantifies the volume of the overlapping region between two fibers and can be normalized in order to provide a similarity index ranging from 1, when fibers are exactly alike, to 0, when they are completely different. Furthermore, we provide the means to calculate the probability Figure 13: Tract probability maps for the averaged white matter fiber bundles over 21 subjects. The bundles have been automatically extracted from each subject individually using the clustering algorithm presented in section 2.3 and the queries presented in table 1. The subjects where previously registered to a white matter fiber atlas and full brain tractography was performed by densely seeding on the whole cerebral white matter. Maximum intensity projection is used for the color intensity. Color code ranges from red, when the probability of the voxel belonging to the bundle is 1.0 to yellow, when the probability of the voxel belonging to the bundle is 0.2. that a point in space belongs to a bundle, the tract probability map of the bundle. Tract probability maps were recently introduced as a tool to perform bundle-specific quantification of scalar quantities like Fractional Anisotropy [START_REF] Hua | Tract probability maps in stereotaxic spaces: Analyses of white matter anatomy and tract-specific quantification[END_REF]. Moreover, tract probability maps have a wide range of applications. For instance, in this work we use them to relate the results of our clustering algorithm with anatomical information given by a white matter atlas. This enables us to perform consistent identification of white matter structures among subjects. These maps could aid as a feature in DTI registration

As part of the clustering, we create a dendrogram, fig. 6 from the full brain tractography. The dendrogram is a versatile tool to represent clustering results, ita a precomputed tree structure where each joint is a candidate cluster. Among other advantages, it allows interactive exploration of the clustering results at different scales and it is resilient to outliers. Moreover, clusters can be selected at different scales using different postprocessing methodologies. In particular, we take advantage of this feature using the tract probability map for Figure 14: Tract probability maps for the averaged white matter fiber bundles over 21 subjects. The bundles have been automatically extracted from each subject individually using the clustering algorithm presented in section 2.3 and the queries presented in table 1. The subjects where previously registered to a white matter fiber atlas and full brain tractography was performed by densely seeding on the whole cerebral white matter. Maximum intensity projection is used for the color intensity. Color code ranges from red, when the probability of the voxel belonging to the bundle is 1.0 to yellow, when the probability of the voxel belonging to the bundle is 0.2.

each cluster and a volumetric atlas with a parcellation of the cerebral gyri. We combine both of these tools with the dendrogram in order to produce a tract query algorithm. For each automatically extracted bundle, the similarity between the mean bundle across all subjects and the bundle extracted from each subject individually was calculated by means of the normalized inner product, eq. 2.2. Similarity ranges from 1, bundles are identical, to 0, bundles are completely different. Boxes span between the second and third quartiles of the similarity with the mean bundle, the red bar is the median similarity. Moreover, the whiskers indicate the bundles whose similarity value with the mean bundle is the smallest and largest within 1.5 times the interquartile distance of the mean similarity and the '+' symbols indicate outliers.

We apply our clustering and querying process to 21 healthy subjects using queries for 13 white matter structures shown in table 1. This produces a set of white matter fiber bundles shown in fig.s 9 to 12. Then, we use our mathematical framework to generate a population-averaged GP and tract probability maps for each identified bundle. These maps allow qualitative analysis. The result of the clustering on the whole population can be analyzed through visual inspection of the population-averaged tract probability maps, fig.s 13 and 14. They allow quantitative analysis: measuring the similarity of the bundle INRIA corresponding to each subject against the population-averaged GP for each bundle, fig. 15, allows us to quantify and analyze the dispersion of the similarity within the population.

Next we look at the qualitative and quantitative assessment of the results of our clustering. Firstly, we note that tracts obtained by means of our clustering-querying procedures (fig.s 9 to 12), are consistent with manually obtained tracts by experts in diffusion MRI images (Wakana et al., 2004, Figures 3, 4, and 5) and macroscopical preparations (Lawes et al., 2008, Figures 5iii and6ii). Since our similarity metric handles partial overlap of fibers, we are able to correctly cluster fibers diverging from complex tracts like the arcuate fasciculus and the cingulum which innervate cortical and subcortical regions. This is observable in fig. 11. Furthermore, visual inspection of the tract probability maps, fig.s 13 and 14, shows that the mean bundles are consistent with probability maps obtained by manual selection of the bundles in Hua et al. (2008, Figure 3) and chemical staining on post-mortem brains [START_REF] Bürgel | White matter fiber tracts of the human brain: Three-dimensional mapping at microscopic resolution, topography and intersubject variability[END_REF]. The high mean similarity with low dispersion on the cingulum, the cortico spinal tract and the inferior fronto-occipital fasciculus is consistent with histological studies [START_REF] Bürgel | White matter fiber tracts of the human brain: Three-dimensional mapping at microscopic resolution, topography and intersubject variability[END_REF]. Additionally, the relatively lower and more disperse mean similarity on the arcuate and uncinate fasciculus is also consistent with previously cited histological studies.

While the method is fully automated, it depends on the parcellation used to produce the tracts. We noticed two cases where non-anatomically coherent low probability regions (yellow) are present: the left section of the fornix, fig. 13(c), and the left uncinate fasciculus, fig. 14(c). In the case of the left side of the fornix, fig. 13(c), even though the high probability region (red) is consistent with the previously cited anatomical studies, there is a low probability region corresponding to the optical nerve that can be seen in this map. Careful analysis of the clustered fibers generating this region shows that this artifact comes from the used tractography techniques; these fibers go all the way from the posterior column of the fornix to the frontal sections of the optical nerve. This is not anatomically correct and is the source of the high variability exhibited by this tract in the quantitative analysis, fig. 15. In the case of the left uncinate fasciculus, fig. 14(c), there is a high probability region (red) which is consistent with the previously cited anatomical studies but there is a low probability region corresponding to a section of the inferior-longitudinal fasciculus. This increases the variability of the clusters corresponding to the left uncinate fasciculus region, fig. 15. Analysis of the clustered fibers constituting this bundle and the query used to obtain the cluster corresponding to the uncinate fasciculus, table 1, shows that this artifact is produced by the white matter parcellation used to create the tracts. The absence of a parcel corresponding to the temporal pole in the atlas leads us to use the medial temporal gyrus for the corresponding query, shown in fig. 7 right side in light green. This allows fibers fully contained within the medial temporal gyrus to be incorporated in the uncinate fasciculus as a result of the corresponding query. In the future, we will investigate the use of alternate atlases or expert-corrected ROIs to remedy this problem.

We presented a mathematical framework to perform statistical operations on white matter fibers obtained from diffusion tensor MRI tractography. We showed that this framework spans an inner product space, and that this inner product space is appropriate for performing statistical operations among white matter bundles. We used this framework to build a clustering algorithm and a tract querying algorithm which allow automatic fiber bundle identification with the white matter queries as a sole parameter. Then, we applied these algorithms to 21 subjects registered to a publicly available atlas with a parcellation of the white matter gyri that we used as a priori anatomical information. Finally, we showed that the results of our clustering were consistent with studies carried by dissection macroscopical preparations and chemical staining. Thus, our framework has proved suitable to perform group studies of the cerebral white matter with minimal user interaction.

Acknowledgments

This work was partly supported by the INRIA ARC Diffusion MRI Program, the CORDI-S INRIA program and the Odyssée-EADS Grant 2118. Ragini Verma and Luke Bloy acknowledge support from the NIH grants R01MH079938 and T32-EB000814 respectively. The data was acquired as part of the NIH grant R01MH060722. Demian Wassermann wants to acknowledge Rebecca Wolpin for English language editing and Romain Veltz for helpful discussions.

A Gaussian Process framework for white matter fibers

Gaussian Processes (GP) are a rich mathematical framework that has been used in a wide variety of research fields. Examples of this being machine learning [START_REF] Mackay | Introduction to gaussian processes[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], spatial statistics [START_REF] Stein | Statistical Interpolation of Spatial Data: Some Theory for Kriging[END_REF] and modeling of observational data [START_REF] Wahba | Spline Models for Observational Data[END_REF]. Furthermore, they are also related to kernel machines and radial basis functions [START_REF] Wahba | Spline Models for Observational Data[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. In this work we take advantage of the GP to produce a parametric representation of fiber bundles. Particularly, we take advantage of the capability to incorporate different types of hypotheses into the model. Moreover, this representation provides a simple way to linearly combine fibers into bundles and to measure similarity through an inner product operation which we develop in section A.3. Moreover, it provides a natural method to calculate tract probability maps, developed in section A.4.

A.1 Single fibers as Gaussian Processes

The parametric representation of a single fiber bundle is based on two hypotheses. Firstly, smoothness: due to the fact that fiber bundles in the brain do not have sharp angles [START_REF] Basser | In vivo fiber tractography using DT-MRI data[END_REF], we consider that the least wiggled trajectory joining the sample points of a fiber INRIA represents that fiber in a most probable manner. Secondly, diffusion associated uncertainty: the uncertainty of the position of each point on the bundle can be modeled using the water diffusion profile at that position [START_REF] Basser | Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI[END_REF].

Using these two hypotheses, we write y F (p), the GP for the indicator function of fiber F, y(p) for clarity, as the combination of two other GPs, z(p) and (p). The process z(p) represents the smoothness of the trajectory in space, its parameters are inferred from the point sequence obtained through the tractography of a fiber: f = {f 1 , . . . , f |f | }⊂R 3 . The process (p) adds a variability to the fiber at p using full diffusion information and it is inferred from the tensor field over the fiber Σ(f 1 ), . . . , Σ(f |f | ). We show now that z(p) and (p) can be modelled as GPs by characterizing them through covariance functions.

A.1.1 Smoothness

We formulate a probabilistic prior for the probability of a function z(•) being the indicator function of a fiber. From the hypotheses that an indicator function with smaller curvature is a more probable representation of a given trajectory and that the original trajectory F is C 2 . We express these hypotheses as the probabilistic prior

-log(P {z(•)}) = 1 2 Ω |D 2 z(p)| 2 dp + const, Ω⊂R 3 (8)
where the linear operator D maps z(•) to its derivative. Then, we characterize z(•) from this prior. To begin with, we rewrite the prior over a finite sampling {p i } i ⊂Ω. Thus, letting [z(p i )] i be the column vector of values of z(•) at each sampled point, we reformulate eq. 8 as

-log(P {[z(p i )] i }) = 1 2 [z(p i )] T i D 2 (D 2 ) T [z(p i )] i + const (9)
which is the p.d.f of a multivariate Gaussian with covariance matrix

S = (D 2 (D 2 ) T ) -1 = (D 4 ) -1 .
Hence, the function z(•) is represented by the GP, z(p)∼GP(z * (p), s(p, p )).

In order to fully characterize the GP for z(p), we need to find an explicit formulation for s(p, p ). We do this by first rewriting S = (D 4 ) -1 as D 4 S = I. This reveals s(p, p ) to be the Green function of the fourth derivative operator [START_REF] Williams | Gaussian process implicit surfaces[END_REF][START_REF] Wahba | Spline Models for Observational Data[END_REF])

Ω D 4 (u, w)s(w, v)dw = δ(u -v). (10)
Next, the solution for eq. 10 inside a sphere in R 3 of radius R is

s(p, p ) := ψ( p -p ),ψ(r) = 2|r| 3 -3Rr 2 + R 3 r≤R 0 r > R , (11)
where the constants have been chosen such that s(p, p ) is a positive semi-definite symmetric function. Thus, we have fully characterized the probabilistic space of functions which describes the blurred indicator function for a smooth trajectory. Furthermore, the covariance function in eq. 11 is equivalent to a bounded thin-plate spline regularizer, which accurately models smooth spatial data [START_REF] Wahba | Spline Models for Observational Data[END_REF].

A.1.2 Diffusion associated uncertainty

We now represent the diffusion information of a point in the trajectory by a blurring GP at that point. The usual practice in GP literature is to model uncertainty at a given point by an isotropic blurring function [START_REF] Mackay | Introduction to gaussian processes[END_REF]. In our case, we use anisotropical blurring at every given point sampled from fiber by means of a second covariance function based on Diffusion Tensor information. This covariance function is built with convolution kernels: Let g(•) be the GP of an isotropic blurring function, then its convolution with a kernel k(•),

(p) = k(w; p)g(w)dw is a GP with covariance function [START_REF] Paciorek | Spatial modelling using a new class of nonstationary covariance functions[END_REF],

u(p, p ) = k(w; p)k(w; p )dw. (12)
We express the blurring kernel k(•; p) as the probability of a particle going from p to w in a time τ in terms of the Diffusion Tensor Σ(p) [START_REF] Basser | MR diffusion tensor spectroscopy and imaging[END_REF]:

k(w; p) = P {w|p, τ, Σ(p)} = = 1 √ (4πτ ) 3 |Σ(p)| e -1 4τ (w-p) T Σ -1 (p)(w-p) .
Then, by performing the integral in eq. 12, we write the covariance function for the anisotropical blurring process as,

u(p, p ) = 1 (4πτ ) 3 |Σ + Σ | e -1 4τ (p-p ) T (Σ+Σ ) -1 (p-p ) ,
where Σ = Σ(p) and Σ = Σ(p ). The blurring at each point is represented by a zero mean GP, (p) ∼ GP (0, u(p, p )) .

A.1.3 Gaussian Process representation of a fiber

Finally, we write the GP formulation for the blurred indicator function of the fiber,

y(p)∼GP y * (p) = z * (p); c(p, p ) = s(p, p ) + u(p, p ) . (13)

INRIA

Up to this point we have a Gaussian Process-based model for the white matter fibers. This model incorporates spatial and diffusion information. Moreover, within this model we can linearly combine fibers into bundles as show in eq. 2 and quantify similarity as shown in eq. 3. In the remainder of this section we show how to characterize the value of the indicator function y(p) at any test point p ∈ R 3 , how to effectively calculate the similarity between two fibers and finally, how to calculate the tract probability map for a bundle.

A.2 Calculating the indicator function value distribution for a test point

We want to calculate a p.d.f. for the value of y(•) at a test point p∈R 3 . This is a simple operation as we have GP representation of the indicator function y(•) for a trajectory F.

due to marginalization properties of the GPs [START_REF] Mackay | Introduction to gaussian processes[END_REF]. Setting y(p) to take the constant value l when p is a point that belongs to the fiber trajectory F, the mean and covariance functions can be calculated in the following way [START_REF] Mackay | Introduction to gaussian processes[END_REF]:

y * (p) = S f (p) T C ff 1l σ 2 (p) = s p -S f (p) T C -1 ff S f (p), (15)
where [S f (p)] i = [s(f i , p)] i , [C ff ] ij = [c(f i , f j )] ij with 1 ≤ i, j ≤ |f |, 1 is the vector with all ones and s p = s(p, p); the functions s(•, •) and c(•, •) were defined in eq. 11 and eq. 13. This formulation is equivalent to "train" a Gaussian Process-based regression with values l at the sampled fiber points and 0 everywhere else. The parameter R is taken as the minimal distance between two consecutive points in f . Examples of blurred indicator functions can be seen in figures 3 and 4.

A.3 Calculating the deterministic inner product between two fiber bundles

In order to produce an inner product space for fiber bundles, we must provide an inner product operation. We have previously introduce this operation in eq. 3. In this section we give a more formal definition and show how to calculate it efficiently. We look for an inner product that quantifies fiber bundle overlapping in a deterministic manner. This simplifies the mathematical treatment of our inner product space, hence its applicability for large scale computations such as clustering of densely sampled full brain tractographies. In order to provide a deterministic inner product, we need to make a decision about the value of y(p) defined in eq. 15. We do this by employing decision theory [START_REF] De Groot | Optimal Statistical Decisions[END_REF]): To make a decision about the value of y(p) at p, we use a point-like prediction, y + (p). is an inner product [START_REF] Schmidt | Über die auflösung linearer gleichungen mit unendlich vielen unbekannten[END_REF]. Furthermore, it can be easily computed by replacing eq. 15 in the previous formula:

F, F := R 3 S f (p) T C ff 1l T S f (p) T C f f 1l dp = R 3 l1 T C T ff S f (p)S f (p) T (C f f 1l) dp.
Moreover, we can simplify the previous expression using that S f (p) and S f (p) are the only vectors depending on p, see section A.2:

F, F := l1 T C T ff R 3 S f (p)S f (p) T dp C f f 1l
where

where h is a bandwidth parameter and σ 2 (p) is defined in eq. 15. Then, the tract probability map for a bundle F on a domain Ω, is calculated by evaluating P {p ∈ F} at every point p ∈ Ω. In order to illustrate the probabilistic map for a bundle, color-coded surfaces and a probability map over an FA image for anatomical bundles are shown in fig. 5. These tract probabilistic maps are highly similar to the hand-obtained ones by [START_REF] Hua | Tract probability maps in stereotaxic spaces: Analyses of white matter anatomy and tract-specific quantification[END_REF] on DTI images and [START_REF] Bürgel | White matter fiber tracts of the human brain: Three-dimensional mapping at microscopic resolution, topography and intersubject variability[END_REF] by means of chemical staining. This two previous works have shown that this tract probability maps are an appropriate tool to perform statistics on white matter fiber bundles.

Figure 1 :

 1 Figure 1: Axons enter and leave anatomical bundles. Image (a) reproduced from O'Donnell (2006) with permission of Jimmy Fallon, UCI. Image (b) manually selected Arcuate fasciculus after full brain tractography.

Figure 2 :

 2 Figure 2: Gaussian Process representation of the fiber: We model each fiber as a blurred indicator function. Each fiber obtained by tractography is constituted as a sample of points over a line, which is shown in red line on the left side of the figure, and its corresponding diffusion tensor field, indicated by the blue tensors over the red fiber. The sampled points are then used to generate a blurred indicator function with a maximal level set which corresponds to the fiber (central image of the figure). Then, the fiber is blurred in accordance with the information provided by the underlying diffusion tensor field. The resulting representation, on the right side of the figure, is a smooth indicator function blurred along the fiber direction.

Figure 3 :Figure 4 :

 34 Figure 3: Mean indicator function for four fiber tracts (a-d) and mean indicator function for the bundle formed by averaging them according to our framework (e). Blue color means that the bundle is more likely to cross that voxel while red color means it is not likely that the bundle traverses that voxel. Fibers were manually selected from a full brain tractography and belong to the Cingulate Cortex section of the Cingulum (CgC). FF & UNC FF & FF FF & FF CgC & CgC

  Figure 5: Iso-probability surfaces for manually selected fiber bundles, fig. (a-c), and tract probability map over FA for an automatically obtained bundle, fig. (d), see section 2.3. The probability at each voxel is calculated using eq. 6. Color code for fig.(a-c) is as follows, Blue: P {p∈B i } = .01, Yellow: P {p∈B i } = .2, Red: P {p∈B i } = .6.

Figure 6 :

 6 Figure 6: Illustration of our clustering algorithm defined in section 2.3. The algorithm starts from all the single-fiber bundles obtained by a full brain tractography, shown at the bottom of the figure, and joins them into multi-fiber bundles according to our similarity measure (eq. 3). This generates a tree structure called dendrogram. Finally every joint is a candidate cluster. Sample clusters are shown on the side of the dendrogram and their positions on the dendrogram are marked with red dots.

Figure 7 :

 7 Figure 7: Illustration of anatomical queries based on volumetric information. A Fractional Anisotropy image and an image of a parcellation of the white matter are shown. The colored stars indicate the anatomical queries. On the left: Purple and red stars tag the inferio-frontal gyrus and middle-occipital gyrus on the left and right sides respectively, this query corresponds to the inferio-fronto-occipital fasciculus on the left and right hemispheres. Yellow stars tag the left and right middle-frontal-orbital gyrus, this query corresponds to the frontal forceps. On the right: white stars tag the pre and post central gyri, the angular, supra-marginal and superio-temporal gyri, this query corresponds to the arcuate fasciculus. The results of this queries on 4 subjects are exhibited in fig.s 9 to 12 and tract probability maps for the mean bundles for every subject are shown in fig.s 13 and 14.

Figure 9 :

 9 Figure9: White matter fiber bundles automatically clustered from a database of 21 subjects using the clustering algorithm presented in section 2.3. The subjects where previously registered to a white matter fiber atlas and full brain tractography was performed by densely seeding on the whole cerebral white matter. Four subjects are shown.

Figure 6 :

 6 Figure 6: Illustration of anatomical queries based on volumetric information. A Fractional Anisotropy image and an image of a parcellation of the white matter are shown. The colored stars indicate the anatomical queries. On the left: Purple and red stars tag the inferio-frontal gyrus and middle-occipital gyrus on the left and right sides respectively, this query corresponds to the inferio-fronto-occipital fasciculus on the left and right hemispheres. Yellow stars tag the left and right middle-frontal-orbital gyrus, this query corresponds to the frontal forceps. On the right: white stars tag the pre and post central gyri, the angular, supra-marginal and superio-temporal gyri, this query corresponds to the arcuate fasciculus. The results of this queries on 4 subjects are exhibited in Figures 7 to 10 and tract probability maps for the mean bundles for every subject are shown in Figures 11 and 12.

Figure 8 :

 8 Figure8: The procedure used to cluster white matter fibers into anatomical bundles and produce the tract probability maps for each bundle: We registered 21 Diffusion Tensor MR images using[START_REF] Yang | Diffusion Tensor Image Registration Using Tensor Geometry and Orientation Features[END_REF]. We performed full brain tractography obtaining around 10.000 fibers per brain. We produced the Gaussian Process representation for each fiber as we describe in section 2.2. We identified anatomical bundles like the arcuate or the uncinate fasciculus by applying our clustering and tract querying algorithms to each subject individually as we describe in section 2.3. Finally we produced a population-averaged Gaussian Process for each identified bundle and the corresponding tract probability map with the methodology described in section 2.2

Figure 9 :

 9 Figure9: White matter fiber bundles automatically clustered from a database of 21 subjects using the clustering algorithm presented in section 2.3 and the queries presented on table 1. The subjects where previously registered to a white matter fiber atlas and full brain tractography was performed by densely seeding on the whole cerebral white matter. Four subjects are shown.

Figure 10 :Figure 11 :

 1011 Figure10: White matter fiber bundles automatically clustered from a database of 21 subjects using the clustering algorithm presented in section 2.3 and the queries presented on table 1. The subjects where previously registered to a white matter fiber atlas and full brain tractography was performed by densely seeding on the whole cerebral white matter. Four subjects are shown.

  Figure12: White matter fiber bundles automatically clustered from a database of 21 subjects using the clustering algorithm presented in section 2.3. The subjects where previously registered to a white matter fiber atlas and full brain tractography was performed by densely seeding on the whole cerebral white matter. Four subjects are shown.

Figure 15 :

 15 Figure15: Quantitative assessment of bundle coherence among subjects. For each automatically extracted bundle, the similarity between the mean bundle across all subjects and the bundle extracted from each subject individually was calculated by means of the normalized inner product, eq. 2.2. Similarity ranges from 1, bundles are identical, to 0, bundles are completely different. Boxes span between the second and third quartiles of the similarity with the mean bundle, the red bar is the median similarity. Moreover, the whiskers indicate the bundles whose similarity value with the mean bundle is the smallest and largest within 1.5 times the interquartile distance of the mean similarity and the '+' symbols indicate outliers.

  The p.d.f. of y(p), given the tractographed point sequence f = {f 1 , . . . , f |f | }⊂R 3 and its corresponding tensor field Σ = {Σ(f 1 ), . . . , Σ(f |f | )}⊂SP D(3) are characterized as Gaussiandistributed random variable. More precisely, y(p) is the Gaussian distribution, (y(p)|f , Σ, p) ∼G y * (p), σ 2 (p)

28DemianF

  Wassermann, Luke Bloy, Efstathios Kanterakis, Ragini Verma and Rachid Deriche This prediction is taken in order to minimise the error in the squared norm induced by our inner product eq. 4,argmin y + (p) y + (p) -y(p) 2 P {y(p)|f , t, p} dy(p) = y * (p). (16)Thus letting the mean value of y(p), y * (p), be an appropriate estimator of the value of y(p) at p. Then, as y * (p) is square integrable due to its definition in eq. 15,

yy

  i , p)s(f j , p)dp can be calculated analytically. Moreover, the similarity among two bundles of white matter fibers, B = {F 1 , . . . , F N } and B = {F 1 , . . . , F M }, F j (p) , can be calculated in a simple manner by using the linearity and symmetry properties of the inner product operation. Then, the previous equation becomes Fi (p), y F j (p) .(17)Examples of the deterministic inner product of two fibers and of the product of two mean indicator functions are shown in fig.4. It can be noted that the inner product quantifies partial fiber overlapping.INRIAA.4 Calculation of the tract probability mapHaving characterized y(•) as a GP in section A.1, we can express the probability that a point p in R 3 is contained in a bundle F, this map is called the tract probability map. We do this by calculating the probability of y(p) = l or, equivalently the concentration of the random value y(p) around l,P {p∈F} := P {y(p) = l|f , Σ, p} ∝E[θ(y(p) -l)|f , t, p].(18)Then, we quantify concentration of y(p) around l by the symmetric kernel θ : R→[0, 1].To ease the equations and the computation time, we take a Gaussian kernel, θ(x), with standard deviation h, θ(y(p) -l) calculate eq. 18 asE[θ(y(p) -l)|f ,t, p] = θ(y(p) -l)P {y(p)|f , t, p} dy(p) which leads to P {p∈F} ∝E[θ(y(p) -l)|f , t, p] = 1 2 π(h 2 + σ 2 (p)) .

Table 1 :

 1 Queries applied to a set of 21 registered subjects in order to automatically extract major white matter fiber bundles using the algorithm presented in section 2.3. The results of this queries on 4 subjects are exhibited in fig.s 9 to 12 and tract probability maps for the tract probability maps corresponding to the population-averaged bundles for every subject are shown in fig.s 13 and 14

	INRIA

The whole process is shown in fig.

8

RR n°7005

RR n°700524Demian Wassermann, Luke Bloy, Efstathios Kanterakis, Ragini Verma and Rachid Deriche


Marcel Hennes

Dominique Demaille

Gilles Patriarche

Thomas Tran

Yunlin Zheng

Franck Vidal
email: franck.vidal@sorbonne-universite.fr

Strain, magnetic anisotropy, and composition modulation in hybrid metal-oxide vertically assembled nanocomposites

Keywords: nanostructure, ferromagnetic, plasmonic, composite, laser ablation

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

Introduction

VANs are nanocomposite systems made of nanowires (NWs) epitaxially coupled to a surrounding matrix (Figure 1a). While early work on these hybrid thin films already demonstrated the possibility to grow embedded metallic nanostructures 1 , most of the subsequent research focused essentially on all oxide VANs (oxide-oxide VANs). Especially the combination of perovskite and spinel oxides has attracted tremendous attention and led to major advances in the development of multiferroic nanocomposites [2][3][4][5][6][7] . In recent years however, the interest for hybrid VANs composed of an oxide matrix and embedded nanowires made of transition or noble metals (metal-oxide VANs) has grown considerably 8- 27 . These systems naturally offer some advantages that may lead to advances in a variety of fields such as nanomagnetism, or plasmonics.

Growth and structural properties of metal-oxide VANs

The growth of epitaxial metal-oxide VANs is based on self-assembly mechanisms driven by metal-oxide phase segregation during synthesis. Pulsed laser deposition is particularly suitable to obtain the desired metal-oxide structures, allowing for deposition at high temperature (>600°C) and under high vacuum conditions (<10 -5 mbar). VANs can be grown from a single target containing oxide components 1 or by using a multi-target approach based on the alternating deposition of the metallic and oxide species 12,14 . Such sequential growth schemes generally rely on the repetition of one basic sequence where submonolayer amounts of the materials are deposited. This offers precise control over the composition. CoNi-SrTiO3 VANs for example can be grown by shooting alternatively on CoO, NiO and SrTiO3 (STO) targets. The NW diameter and volume fraction in the epilayer can then be controlled by the ratio of the number 3 of shots on CoO-NiO targets and on the STO target 12,14 . The Co content in the NWs is tuned by the ratio of the shots on CoO and NiO. The epilayer thickness, thus the NW length, results from the total number of basic sequences. With this procedure, arrays of ultrathin (diameter < 5 nm) Co, Ni and CoNi alloy NWs with high densities (> 10 11 cm -2 ) have been obtained. VAN growth mechanisms are complex, as they involve simultaneous phase segregation, nucleation and coalescence phenomena in oxide and metallic phases. However, the initial stages of VAN self-assembly can be grasped by using appropriate computer simulation approaches. Ni-STO VAN thin film deposition was for example successfully described using a coarse-grained kinetic Monte-Carlo scheme 21 which allowed to identify the parameters affecting the final assembly. While low growth temperatures result in short, isolated Ni nanoinclusions in the matrix, increasing the deposition temperature beyond 820 K favors the growth of rods. Systematic increase of the Ni/STO ratio in a basic sequence yields larger average NW diameters, while the overall wire density remains constant. The relaxation time between deposition steps also constitutes a key parameter to accurately tune the NW size and density.

One central feature of VANs is the epitaxial NW/matrix coupling along the vertical interface resulting in strained NWs 28 , which permits to overcome the conventional critical thickness for strain relaxation in planar heterostructures, as discussed by MacManus-Driscoll et al. and Chen et al. in oxide-oxide VANs [29][30] .

In metal-oxide VANs, the dilation along the nanowire backbone can reach several percent when choosing appropriate metal/oxide combinations: In Ni-CeO2 VAN thin films 13 for example, an axial tensile strain of 1.5% was reported. In highly lattice-mismatched systems such as CoNi-(BaSr)TiO3 14 or CoNi-STO 27 , the latter can even reach up to 4%, as shown in Figure 2a and 2b, while the radial compressive strain remains negligible.

As the NW/matrix interface energy scales with the NW diameter (~d) while the volume elastic energy scales with its sectional area (~d 2 ), the axial strain in NWs increases when the diameter decreases (Figure 2a). Calculations using a Frenkel-Kontorova toy model 13 indicate that the tensile strain state observed in NWs is metastable and that the sequential growth scheme, with progressive adaptation of NW atoms to local variations of the periodic matrix potential, plays an important role in stabilizing the large axial strains. However, the actual strain state is more complex than what can be predicted using this simplified 1D approach which only takes the horizontal (001) lattice planes into account. In fact, (101) and (111) dense planes also play an important role in imposing a lattice and orientation matching between segments of a NW and the matrix 22,27 . Beyond the mean axial strain observed, partially relaxed segments were identified in NWs, with (101) and (111) planes almost parallel to those of the matrix and further relaxed segments in Ni NWs, with a small rotation around [100]-type axis. Figure 2c shows a cross-sectional transmission electron microscopy image of a Ni-STO VAN, unraveling local tilts of 4°. Such a rotation leads to lattice matching along the NW/matrix vertical heterointerface for one set of (101) planes, but it also results in an increase of structural disorder 22 . In ultrathin wires, these complex structural rearrangements ultimately lead to an almost complete loss of crystalline structure 22 . Strain relaxation mechanisms in metal/oxide VANs thus differ from their all oxide counterparts, where substitutional defect formation, misfit dislocations and antiphase boundaries were often reported 31 .

Tailoring the magnetic anisotropy in metal-oxide VANs

The magnetic anisotropy K is a key property of nanomagnets because it determines the direction of easy magnetization and the temperature range of stability of the magnetization. VANs are ideal systems for tailoring the magnetic anisotropy. Indeed, as shown schematically in Figure 1, by tuning the growth parameters, and thus the final nanoarchitecture of the VANs, it is possible to control the distinct sources of K independently.

The shape of the embedded nanopillars offers direct control over the magnetostatic part of K, the so-called shape anisotropy Kms. This contribution scales with the square of the magnetization Ms and is given, in the case of nanowires with large aspect ratio, by: Kms =µ0Ms 2 /4.

(1)

Thus, using materials with a large value of Ms is a simple way to obtain a large uniaxial anisotropy. This has been shown for Fe nanopillars obtained by decomposition, as reported in the seminal work of Ref. 1. In such metal-oxide VANs, nanowires exhibited a robust ferromagnetic response at room temperature with a clear easy axis along their backbone and a large coercive field for aspect ratios larger than 5. More recently, the possibility to tune the shape of pillars, in particular their cross section, by adjusting the growth rate, was demonstrated for Fe-BaTiO3 VANs 32 resulting in a tunable magnetization and coercivity.

The epitaxial relationship determines the orientation of the embedded nanocrystal and thus, the easy magnetic directions favored by the magnetocrystalline anisotropy Kmc. For most metals studied so far, this contribution has only played a minor role 27 and it remains to be fully elucidated, whether more complex metallic nanoalloys such as CoPt and FePt, exhibiting hard magnetic phases as a result of chemical short-range ordering, can be integrated into VANs.

Finally, the large axial strain within the nanopillars is a lever that can be used to tune the magnetoelastic anisotropy Kme. For a cubic structure and [001] direction parallel to the longitudinal axis, it is given by:

Kme =3λ001 (c11-c12)(εzz-εrr)/2, ( 2
)
λ001 is the magnetostriction coefficient, cij are elastic constants, εzz and εrr are the axial and radial strain in the nanowires, respectively. In the case of VANs containing alloy nanowires (such as CoxNi1-x), one benefits from an additional degree of freedom: the composition of the metallic inclusions, adjusted precisely via the growth sequence: In pure Ni, λ001 is negative. In CoxNi1-x alloys, λ001 increases with x, reaches 0 around x=18% and is positive above this composition.

Thus, by controlling both, the strain and the composition of CoxNi1-x nanowires, it becomes possible to obtain à-la-carte uniaxial magnetic anisotropies 12,14,24 .

In pure Ni nanowires, there is a threshold value, εcrit, for the axial strain that determines the transition from an easy magnetic axis along the axis of the wires to a hard axis. The value of εcrit is determined by the competition of the shape anisotropy (equation 1) and the magneto-elastic anisotropy (equation 2).

For Ni nanowires embedded epitaxially in a CeO2 matrix (grown on SrTiO3(001) substrate), the axial dilation is of the order of 1%, near εcrit, and the system is nearly isotropic from a magnetic point of view 12 . In a SrTiO3 matrix (grown on SrTiO3(001) substrate), the strain in Ni nanowires is much larger, as shown in Figure 2a and one can expect a reversal of the magnetic behavior: the axis of the wires switches from easy to hard. This is indeed what has been observed 24 , as demonstrated in Figure 3a-b. In Ni-STO VANs, the strain largely exceeds εcrit and the axis of the wires becomes the hard axis with a minimal value of the remanence, as shown in Figure 3b.

In VANs incorporating CoxNi1-x nanowires, the axial strain is still large (see Figure 2b). Based on the above considerations, for x close to 18%, only the shape anisotropy is present in the system. This could indeed be verified by quantitative measurements of the magnetic anisotropy in Co0.2Ni0.8-(BaSr)TiO3

VANs through a combination of cyclic magnetometry (figure 3c) and ferromagnetic resonance measurements 14 . For higher Co contents, the magnetoelastic contribution should lead to an enhancement of the uniaxial anisotropy. This was confirmed in CoxNi1-x-CeO2 VANs: Figure 3d shows magnetic hysteresis cycles recorded at room temperature for a Co0.63Ni0.37-CeO2

VAN exhibiting a ferromagnetic response with strongly anisotropic behavior.

Extraction of the values of the effective uniaxial anisotropy, Keff, of the system shows that it indeed exceeds the shape anisotropy contribution at high Co content (Figure 3e).

The strain engineering approach thus enables to boost K and to push the superparamagnetic limit far beyond room temperature, even for NWs with diameters smaller than 5 nm. This is an exciting perspective for future data storage applications and might be a key step towards controlling anisotropic transport properties in VAN-based nano-devices 18 .

Towards axial structuration and hybrid nanowires

Metal-oxide VANs with two distinct embedded metallic phases can yield additional appealing functionalities. Figure 4a illustrates some hypothetical nanoarchitectures, reminiscent of the multiple configurations adopted by bimetallic nanoparticles 33,[START_REF] Ferrando | Structure and Properties of Nanoalloys[END_REF] . To what extent can these structures be predicted in real-world VANs? Simulation approaches have been developed during the last two decades and turned into powerful tools used to analyze bi-and trimetallic particle structures [START_REF] Ferrando | [END_REF] . However, the integration of metal-oxide interactions into the calculations 36 remains a great challenge. So far, the theoretical analysis of phase separation tendencies in bimetallic VANs has remained limited to basic thermodynamic considerations 23 .

Embedded structures with core-shell (CS) or multishell architecture provide an in-plane modulation of the wire properties. In freestanding NWs, such configurations have been extensively studied and result from sequential growth approaches, where the core of the NW is obtained first, then covered by the shell material. This is hardly compatible with VAN growth, where the wire and matrix assemble simultaneously. The growth of core-shell VANs is thus challenging.

Most alternative approaches developed so far, such as post-growth annealing, yielded limited success. Partial interdiffusion at the NW/matrix vertical interface has been observed in oxide-oxide VANs 37 giving rise to a "diffuse" CS configuration. Spontaneous decomposition and segregation processes, occurring during growth, have been shown to lead to CS configurations (Cu core-SrO shell 38 and TiO core-MgO shell 39 ) but, to the best of our knowledge, this has not been achieved with two embedded metallic phases.

With regard to their magnetic and charge/spin transport properties, backbone-concentration modulated magnetic NWs appear even more appealing than their CS counterparts. Combining for example soft and hard magnetic materials would allow the growth of self-organized spin-valves with unprecedented small diameters. At first sight, these structures seem less challenging to synthesize. One might for example expect to obtain multilayered NWs (Figure 4a) upon sequential deposition of different metallic species.

Unfortunately, the rather high growth temperatures required to achieve metal/oxide phase separation concomitantly induce pronounced vertical atom mobility. Working with miscible metals requires finding a compromise between in-plane 40 and out-of-plane diffusion 21 , in order to avoid creation of a bimetallic solid solution. Mixing of the components can be avoided by choosing appropriate material combinations. Au-Co or Ag-Co for example remain fully immiscible in the solid phase, even in nanoscale structures. However, the equilibrium VAN configurations resulting from co-deposition of these metals are also strongly influenced by differences in their surface energy. Co-growth of Au-Co-STO and Ag-Co-TiO2 produced bilayered NWs, with the component with lower surface energy located at the NW top 23,41,42 : Au and Ag act as metallic surfactants, pushing incoming Co atoms to the bottom of the NW (Figure 4b), which additionally minimizes the interface between the two immiscible metallic phases.

Similar processes have also been observed in three-phase metal/oxide/oxide composites 43 .

While this impedes the growth of multilayer NWs, the aforementioned noble metals -ferromagnetic bilayer configurations can still be considered highly promising for a variety of applications, especially with regard to their plasmonic functionalities. The potential of self-assembly to yield embedded metallic columns with tunable optical resonances has been recognized early on 44,45 and was recently renewed for VANs 16,46 . Coupling of these plasmonic features with magnetic properties might turn VANs into a beautiful playground for magnetophotonic studies at the nanoscale 47 , especially in an ultrafast regime 48 , where fascinating non-equilibrium phenomena such as sub-ps demagnetization 49 and all optical switching 50 can be triggered by optical fs-pulses.

Summary and outlook

Metal-oxide VANs constitute a novel class of hybrid thin films of special interest in the field of nanomagnetism, providing the opportunity to obtain à-la-carte magnetic anisotropies. Among the challenges ahead, the successful growth of

Figure 1 .

 1 Figure 1. Schematic illustration of a VAN and of the possibilities of magnetic anisotropy control.

Figure 2 .

 2 Figure 2. (a) Axial strain <εzz> as a function of NW diameter in Ni-CeO2 VANs (disks, adapted from Ref. 13 by permission from Springer Nature Customer Service Centre GmbH: Springer Nature, Nano Research, Huge metastable axial strain in ultrathin heteroepitaxial vertically aligned nanowires, V. Schuler, F. J. Bonilla, D. Demaille, A. Coati, A. Vlad, Y. Garreau, M. Sauvage-Simkin, A. Novikova, E. Fonda, S. Hidki, V. Etgens, F. Vidal, Y. Zheng, © 2015) and Ni-STO VANs (triangles, adapted from Ref. 22 with permission). (b) Out-of-plane lattice parameters az of CoxNi1-x-STO VANs, with the cyan line as guide to the eye and the red line representing the bulk CoNi alloy reference [Ref. 27] -Reproduced by permission of The Royal Society of Chemistry. (c) Cross section transmission electron microscopy image of Ni nanowires in SrTiO3 matrix; insets: Fast Fourier Transform of selected areas delimited by light squares. The inset at the top evidences the local misorientation α of the Ni lattice with respect to the SrTiO3 lattice (reproduced from Ref. 22 with permission).

Figure 3 .

 3 Figure 3. (a-c) Magnetic hysteresis cycles at 10 K in Ni-STO (a), in Ni-CeO2 (b) and in Co0.2Ni0.8-(BaSr)TiO3 (c). (d) Magnetic hysteresis cycles at 300 K in a Co0.63Ni0.37-CeO2 VAN. (e) Black line: Magnetostatic anisotropy at low temperature as a function of the Co content in CoxNi1-x alloy. Black triangle: magnetic anisotropy of Ni wires in Ni-CeO2 VAN. Keff, magnetic anisotropy extracted from magnetometry, of Co0.32Ni0.68 (green plain square, T = 10 K; open square, T = 300 K) and Co0.63Ni0.37 (dark red disk, T = 10 K; circle, T = 300 K) nanowires assemblies embedded in CeO2. (a) adapted from Ref. 24 with

Figure 4 .

 4 Figure 4. (a) VANs with nanowires exhibiting core-shell (CS) multilayer, and bilayer segregation patterns. (b) Hybrid Au-Co-SrTiO3 VAN: STEM-HAADF and EDX chemical mapping of Co and Au in cross-sectional view. Reproduced from Ref. 23 with permission.

modulated multi-metallic VANs structures, combining distinct physical properties, will be of utmost interest in the emerging field of magnetoplasmonics. Bimetallic VANs would be an ideal playground to study novel synergy effects on ultrashort timescales, thereby paving the way for an improved photon-based control of magnetic properties in ferromagnetic nanomaterials.