Bertrand Thirion
email: bertrand.thirion@inria.fr

Ga Ël Varoquaux

Elvis Dohmatob

Jean-Baptiste Poline

Which fMRI clustering gives good brain parcellations?

Keywords: Functional neuroimaging, Brain atlas, clustering, Model selection, cross-validation, group studies

Analysis and interpretation of neuroimaging data often require one to divide the brain into a number of regions, or parcels, with homogeneous characteristics, be these regions defined in the brain volume or on on the cortical surface. While predefined brain atlases do not adapt to the signal in the individual subjects images, parcellation approaches use brain activity (e.g. found in some functional contrasts of interest) and clustering techniques to define regions with some degree of signal homogeneity. In this work, we address the question of which clustering technique is appropriate and how to optimize the corresponding model. We use two principled criteria: goodness of fit (accuracy), and reproducibility of the parcellation across bootstrap samples. We study these criteria on both simulated and two task-based functional Magnetic Resonance Imaging datasets for the Ward, spectral and K-means clustering algorithms. We show that in general Ward's clustering performs better than alternative methods with regard to reproducibility and accuracy and that the two criteria diverge regarding the preferred models (reproducibility leading to more conservative solutions), thus deferring the practical decision to a higher level alternative, namely the choice of a trade-off between accuracy and stability.

INTRODUCTION

Brain parcellations divide the brain's spatial domain into a set of non-overlapping regions or modules that show some homogeneity with respect to information provided by one or several image modalities, such as cyto-architecture, anatomical connectivity, functional connectivity, or task-related activation. Brain parcellations are therefore often derived from specific clustering algorithms applied to brain images. Such approaches are generally useful because the voxel sampling grid of the reference space, e.g. the MNI template, is most often at a higher resolution than the brain structures of interest, or at a scale that is too fine for the problem under investigation, yielding an excessive number of brain locations and correlated data. In other words, the structures of interest are rarely at the level of a specific voxel, but at the level of many voxels constituting a (possibly small) brain region. Three strategies are commonly used to study function beyond the voxel description: i) the use of anatomical or functional regions of interest (ROIs), ii) the use of a brain atlas, or iii) the use of data-driven parcellations.

ROI-based analysis has been advocated as a way to focus data analysis on some structures of interest and consists in building a summary of the signal in a predefined region [START_REF] Nieto-Castanon | Region of interest based analysis of functional imaging data[END_REF]. The choice of the region(s) can be based on prior experiments (e.g. [START_REF] Saxe | Divide and conquer: a defense of functional localizers[END_REF]). Note that in extreme cases, the region can reduce to a single voxel, one reported in previous literature as the peak coordinate of a contrast image 1 . The obvious limitation of ROI-based analysis is that the signal present outside the region under consideration is ignored a priori; as a consequence, the results depend heavily on the choice of this ROI, which may not fit well the new data. In the hypothesis testing framework, the smaller number of tests performed may however increase the power of the analysis.

Brain atlases come into play to provide a set of ROIs that cover the brain volume (among many others see e.g. [START_REF] Mazziotta | A probabilistic atlas and reference system for the human brain: International consortium for brain mapping (icbm)[END_REF]; [START_REF] Tzourio-Mazoyer | Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the mni mri single-subject brain[END_REF]; [START_REF] Shattuck | Construction of a 3d probabilistic atlas of human cortical structures[END_REF]). An atlas generally accounts for a certain state of the knowledge of the brain structures (anatomically, functionally or based on connectivity), from which well-defined entities can be distinguished. In other words, an atlas represents a certain labeling of brain structures. Often this labeling is linked to an ontology representing the current knowledge [START_REF] Cieslik | Is there "one" dlpfc in cognitive action control? evidence for heterogeneity from co-activation-based parcellation[END_REF][START_REF] Eickhoff | Co-activation patterns distinguish cortical modules, their connectivity and functional differentiation[END_REF]. In spite of their obvious usefulness, existing atlases are limited in two regards: i) There exist currently many different atlases, but they are mutually inconsistent [START_REF] Bohland | The brain atlas concordance problem: quantitative comparison of anatomical parcellations[END_REF]; ii) A given atlas may not fit the data well. Atlas misfits can be due to image characteristics and processing strategies that have evolved since the atlas creation, or because a given study deals with a population that is not well represented by the subjects used to construct the atlas, or because the information of interest is simply not mapped properly in the given atlas. Atlas misfit is often pronounced with regards to mapping brain function; for instance most anatomical atlases have large frontal brain regions that many researchers would rather divide into smaller ones with more precise functional roles.

Unlike brain atlases, also used to define regions of interest, brain parcellations are data-driven. They do not reflect a pre-defined ontology of brain structures -known anatomical names and concepts-, but they may much better represent the measurements or features of interest, i.e. they provide a better model of the signal [START_REF] Flandin | Improved detection sensitivity in functional MRI data using a brain parcelling technique[END_REF]Simon et al., 2004;[START_REF] Thirion | Dealing with the shortcomings of spatial normalization: multi-subject parcellation of fmri datasets[END_REF][START_REF] Lashkari | Discovering structure in the space of fMRI selectivity profiles[END_REF][START_REF] Lashkari | Search for patterns of functional specificity in the brain: a nonparametric hierarchical bayesian model for group fMRI data[END_REF]. The (anatomical) labeling of these parcels can then be performed with the most appropriate atlas.

While functional parcellations can be used in different contexts, we focus here on finding a well-suited model to obtain local averages of the signal for group studies. These parcel averages can be thought of as a data reduction adapted to various tasks, such as the estimation of brain-level connectivity models (see e.g. [START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF]; [START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF]), of physiological parameters [START_REF] Chaari | Hemodynamic-informed parcellation of fmri data in a joint detection estimation framework[END_REF], for group analysis [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: multi-subject parcellation of fmri datasets[END_REF], the comparison of multiple modalities [START_REF] Eickhoff | Co-activation patterns distinguish cortical modules, their connectivity and functional differentiation[END_REF] or in multivariate models [START_REF] Michel | A supervised clustering approach for fMRI-based inference of brain states[END_REF]. This is especially useful for the analysis of large cohorts of subjects, because this step can reduce the data dimensionality by several orders of magnitude while retaining most of the information of interest. We will show in this paper that common brain atlases, merely reflecting sulco-gyral anatomy, are not detailed enough to yield adequate models of the (functional) data.

Data-driven parcellations can be derived from various image modalities reflecting different neurobiological information, for instance T1 images with anatomical information, such as gyrosulcal anatomy [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on mri scans into gyral based regions of interest[END_REF][START_REF] Klein | 101 labeled brain images and a consistent human cortical labeling protocol[END_REF], post-mortem in vitro receptor autoradiography for cyto-architecture [START_REF] Eickhoff | Organizational principles of human visual cortex revealed by receptor mapping[END_REF][START_REF] Fischl | Cortical folding patterns and predicting cytoarchitecture[END_REF], anatomical connectivity [START_REF] Roca | Inter-subject connectivity-based parcellation of a patch of cerebral cortex[END_REF] with diffusion imaging, or functional features with BOLD data. In this work, we focus on the latter, that we call functional parcellations. These parcellations are currently derived either from resting-state functional Magnetic Resonance Images (rs-fMRIs) [START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF][START_REF] Blumensath | Resting-state FMRI single subject cortical parcellation based on region growing[END_REF][START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF][START_REF] Kahnt | Connectivity-based parcellation of the human orbitofrontal cortex[END_REF][START_REF] Wig | Parcellating an individual subject's cortical and subcortical brain structures using snowball sampling of resting-state correlations[END_REF], from activation data [START_REF] Flandin | Improved detection sensitivity in functional MRI data using a brain parcelling technique[END_REF][START_REF] Lashkari | Discovering structure in the space of fMRI selectivity profiles[END_REF][START_REF] Lashkari | Search for patterns of functional specificity in the brain: a nonparametric hierarchical bayesian model for group fMRI data[END_REF][START_REF] Michel | A supervised clustering approach for fMRI-based inference of brain states[END_REF], or from meta-analyses [START_REF] Eickhoff | Co-activation patterns distinguish cortical modules, their connectivity and functional differentiation[END_REF]. To investigate which parcellations are most appropriate, we restrict our work to activation data that 1 Often, a small sphere will be drawn around this position to average signals locally. This is a provisional file, not the final typeset article have more tractable probabilistic models than resting-state data. We also omit edge-based parcellation methods, such as those described in [START_REF] Wig | Parcellating an individual subject's cortical and subcortical brain structures using snowball sampling of resting-state correlations[END_REF] or [START_REF] Cohen | Defining functional areas in individual human brains using resting functional connectivity mri[END_REF]: while these are certainly useful to segment the cortical surface by revealing abrupt changes in the functional connectivity patterns when crossing region boundaries, they do not lend themselves to model selection due to the absence of a probabilistic framework. This family of approaches is certainly an interesting competitor for future analyzes of functional parcellations performed on the cortical surface.

The most popular parcellation techniques are mixture models [START_REF] Golland | Detection of spatial activation patterns as unsupervised segmentation of fmri data[END_REF][START_REF] Lashkari | Discovering structure in the space of fMRI selectivity profiles[END_REF][START_REF] Lashkari | Search for patterns of functional specificity in the brain: a nonparametric hierarchical bayesian model for group fMRI data[END_REF][START_REF] Tucholka | Probabilistic anatomo-functional parcellation of the cortex: how many regions?[END_REF], variants of the k-means algorithm [START_REF] Flandin | Improved detection sensitivity in functional MRI data using a brain parcelling technique[END_REF][START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF][START_REF] Kahnt | Connectivity-based parcellation of the human orbitofrontal cortex[END_REF], hierarchical clustering [START_REF] Eickhoff | Co-activation patterns distinguish cortical modules, their connectivity and functional differentiation[END_REF][START_REF] Michel | A supervised clustering approach for fMRI-based inference of brain states[END_REF][START_REF] Orban | The richness of taskevoked hemodynamic responses defines a pseudohierarchy of functionally meaningful brain networks[END_REF] and variants thereof [START_REF] Blumensath | Resting-state FMRI single subject cortical parcellation based on region growing[END_REF], spectral clustering [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: multi-subject parcellation of fmri datasets[END_REF][START_REF] Chen | Inferring group-wise consistent multimodal brain networks via multi-view spectral clustering[END_REF][START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF] and dense clustering [START_REF] Hanson | Dense mode clustering in brain maps[END_REF]. Some of these approaches, but not all, impose spatial constraints on the model, and therefore provide spatially connected spatial components. In the multi-subject setting, some models adapt the spatial configuration to each subject (e.g. [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: multi-subject parcellation of fmri datasets[END_REF]; [START_REF] Lashkari | Discovering structure in the space of fMRI selectivity profiles[END_REF][START_REF] Lashkari | Search for patterns of functional specificity in the brain: a nonparametric hierarchical bayesian model for group fMRI data[END_REF]), but most approaches do not.

Parcellations can also be obtained from dictionary learning techniques such as independent components analysis (ICA) and variants of principal components analysis (PCA) [START_REF] Kiviniemi | Functional segmentation of the brain cortex using high model order group pica[END_REF][START_REF] Varoquaux | A group model for stable multi-subject ica on fmri datasets[END_REF][START_REF] Varoquaux | Multi-subject dictionary learning to segment an atlas of brain spontaneous activity[END_REF][START_REF] Varoquaux | Cohort-level brain mapping: learning cognitive atoms to single out specialized regions[END_REF][START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF]. These rely on a linear mixing approach that changes the nature of the problem and implies other probabilistic models.

While parcellation techniques have great potential and can serve as the basis of many further analyses, it is important to assess their relative performance. To the best of our knowledge, no systematic comparison of parcellation methods has been carried out in previous work.

The comparison between clustering techniques is only relevant if for each technique the best possible model is selected. It turns out that model selection for clustering is a notoriously difficult problem, as is any unsupervised problem in which one wishes to identify some structure in noisy data. While in practice the choice of the model may depend on the context of the study (for instance, fitting a given target of interest using region-based signal averages [START_REF] Ghosh | Predicting treatment response from resting state fmri data: Comparison of parcellation approaches[END_REF]), here we derive general rules to compare parcellation models from empirical observations. In the context of brain mapping, two criteria are particularly relevant for model selection: i) the goodness of fit or accuracy of a model, i.e. the ability of the parcellation extracted to model properly the signals of interest on observed and unobserved data, and ii) stability, i.e. the consistency of the parcellations obtained from different sub-groups of a homogeneous population. Importantly, there is a priori no reason why these two criteria should give consistent answers. There have been few attempts to tackle this, such as [START_REF] Tucholka | Probabilistic anatomo-functional parcellation of the cortex: how many regions?[END_REF][START_REF] Kahnt | Connectivity-based parcellation of the human orbitofrontal cortex[END_REF]; [START_REF] Ghosh | Predicting treatment response from resting state fmri data: Comparison of parcellation approaches[END_REF], but these approaches did not model the multi-subject nature of the signal; moreover [START_REF] Tucholka | Probabilistic anatomo-functional parcellation of the cortex: how many regions?[END_REF] were subdividing prior gyrus definition (hence not brain-wide) and they did not benchmark different clustering techniques. In the present work, we present experiments on simulated and real data using different clustering techniques and proper accuracy and reproducibility criteria. To make this tractable computationally and to obtain clear interpretation, we limit ourselves to nonlinear mixing models, i.e. clustering approaches. Note that methods comparison for clustering versus linear mixing models (ICA, variants of sparse PCA) has been addressed e.g. in [START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF], while model order selection for linear model-based region extraction is still an open problem. For similar reasons, we consider the case in which parcels are identical for all subjects.

The remainder of the paper is organized as follows: in Section 2, we introduce the methods tested in this work and the criteria for model evaluation; in Section 3 we describe our experiments on simulated and real data, the results of which are given in Section 4. Conclusions on the choice of optimal processing algorithms and the selection of parcellation schemes are drawn in Section 5.

MATERIAL & METHODS

NOTATION

We start with a given set of n functional images that represent e.g. different contrasts in a given group of subjects. We denote N to be the number of subjects and F the number of functional images (here contrasts) per subject, such that n = N F . These images are typically the results of first-level analysis (standardized effects) and are sampled on a grid of Q voxels. Starting from n fMRI volumes Y = [y 1 , . . . , y Q ] ∈ R n×Q that consist of Q voxels, we seek to cluster these voxels so as to produce a reduced representation of Y.

CLUSTERING METHODS FOR BRAIN FUNCTIONAL PARCELLATION

K-means Algorithm K-means is arguably the most used clustering technique for vector data. It consists of an alternate optimization of i) the assignment u k-means of samples to cluster and ii) the estimation of the cluster centroids.

∀j ∈ [1, Q], u k-means (j) = argmin c∈[1,...,K] Y c -y j
(1)

Y c 1 |c| u k-means (j)=c y j (2)
It explicitly minimizes the inertia, i.e. the sum of squared differences between the samples and their representative cluster centroid. We introduce an approximation for the sake of efficiency: the whole set of feature data used in clustering (several contrasts from all the subjects) of dimension n = N (subjects) × F (contrasts) is reduced by PCA to m = 100 components prior to clustering, capturing about 50% of the variance. It is important to note that k-means clustering of fMRI data are used without explicitly considering their spatial structure, although spatial smoothing prior to clustering can indirectly provide spatial regularization.

Ward's Algorithm As an alternative, we consider a hierarchical agglomerative clustering [START_REF] Johnson | Hierarchical clustering schemes[END_REF]. These procedures start with every voxels x j representing singleton clusters {j} and, at each iteration, a pair of clusters, selected according to a criterion discussed below, is merged into a single cluster. This procedure yields a hierarchy of clusters represented as a binary tree T , also often called a dendrogram [START_REF] Johnson | Hierarchical clustering schemes[END_REF], where each non-terminal node is associated with the cluster obtained by merging its two children clusters.

Among different hierarchical agglomerative clustering procedures, we use the variance-minimizing approach of Ward's algorithm [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF]. In short, two clusters are merged if the resulting cluster minimizes the sum of squared differences of the fMRI signal within all clusters. More formally, at each step of the procedure, we merge the clusters c 1 and c 2 that minimize

∆(c 1 , c 2 ) = j∈c 1 ∪c 2 y j -Y c 1 ∪c 2 2 2 - j∈c 1 y j -Y c 1 2 2 + k∈c 2 y k -Y c 2 2 2 = |c 1 ||c 2 | |c 1 | + |c 2 | Y c 1 -Y c 2 2 2 , (3)
where Y c is the average vector defined in (2). In order to take into account the spatial information, we also add connectivity constraints in the hierarchical clustering algorithm, so that only neighboring clusters can be merged together. In other words, we try to minimize the criterion ∆(c 1 , c 2 ) only for pairs of clusters that share neighboring voxels. Given a number of parcels K, we stop the construction of the tree at the (Q -K) th iteration and retain the corresponding assignment u ward . Note that the data are subject to the same PCA procedure as for k-means clustering.

Spectral clustering Spectral clustering [START_REF] Shi | Automatized clustering and functional geometry of human parietofrontal networks for language, space, and number[END_REF][START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] consists in performing k-means clustering on a representation of the data that preserves the spatial structure yet represents the functional features' similarity2 . This representation is typically obtained by using the first eigenvectors of the Laplacian matrix of the graph that encodes the spatial relationships weighted by the functional features similarity between adjacent locations. For all voxel pairs (i,

j) ∈ [1 • • • Q] 2 , Let W ij = exp(-y i -y j ) 2 2σ 2 f
) if i and j are neighbors

0 otherwise (4)
where we used σ 2 f = mean i∼j y iy j 2 , where the averaging is performed over all pairs of adjacent voxels. σ 2 f is thus the average squared distance between the data across neighboring voxels. W is therefore an adjacency matrix weighted by the functional distance between voxels. We denote ∆ W the diagonal matrix that contains the sum of the rows of W .

Then, let (ξ 1 , .., ξ m ) the first m solutions of W ξ = λ∆ W ξ. The spectral clustering of the dataset is defined as:

u spectral = k-means([ξ 1 , .., ξ m ]), (5)
m = 100 in our experiments. We also tried different (larger or smaller) values, but those did not yield significantly better solutions.

Geometric clustering

To provide a reference for comparison, we also use a clustering algorithm that does not take into account the functional data, but only the spatial coordinates of the voxels. In practice, it is obtained through a k-means clustering of the spatial coordinates.

A MIXED-EFFECTS MODEL OF THE SIGNAL WITHIN PARCELS

We introduce a probabilistic model of the signal of the voxels in a given (fixed

) parcel P k , k ∈ [1, • • • , K],
that includes a random subject effect. Let us first assume that we work with one functional image (F = 1).

Let p be the number of voxels in the parcel, pooled across subjects: it is the size of the parcel multiplied by N ; let y be a p-dimensional vector that denotes the scalar signal in the voxels contained in P k , concatenated across subjects; we model it though the following mixed-effects model:

y = µ1 + Xβ + ε, ( 6
)
where µ is the average signal within the parcel, 1 is a vector of ones of length p, β is a vector of subjectspecific random effects parameters, X the (known) matrix that maps subjects to voxels: for each row, a one in the i th column indicates that the value is from subject i. ε represents the intra-subject variability of the signal within a parcel. It is further assumed that ε and β are independent, normal and centered at 0, with variance σ 2 1 and σ 2 2 that express respectively the within and between subject variance. The probabilistic model of y is thus:

y ∼ N (µ1, σ 2 1 I + σ 2 2 XX T ), (7)
where I is the p × p identity matrix.

The generalization to non-scalar images (for instance, F > 1 images per subject) is obtained by assuming the independence of the observations conditional to the parcellation, hence it decouples into multiple (F ) scalar models. The estimation of the parameters (µ, σ 1 , σ 2 ) is carried out in each parcel

P k , k ∈ [1, • • • , K]
using the maximum likelihood principle; we use an Expectation-Maximization algorithm to estimate the model parameters [START_REF] Meng | Fast em-type implementations for mixed effects models[END_REF].

MODEL SELECTION FOR FUNCTIONAL PARCELLATIONS

A problem that comes naturally with clustering algorithms is the choice of the number K of clusters to be used in the model. To guide this choice we consider four standard measures: BIC, cross-validated likelihood, adjusted rand index, and normalized mutual information.

Bayesian Information Criterion, BIC The goodness of fit of a probabilistic model is given by the loglikelihood of the data and the quality of the model is easily measured using the BIC criterion [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF], that penalizes the negative log-likelihood by the number of parameters used. Within a given parcel P k , this yields the following:

bic(k) = -2 log N (y; µ1, σ 2 1 I + σ 2 2 XX T ) + 3 log (p), ( 8
)
Where 3 is the number of parameters of the model (µ, σ 1 , σ 2 ). Note that all the quantities in this formula (y, µ, σ 2 1 , σ 2 2 , X, p) depend on k, the index of the parcel. bic(k) is summed across parcels in order to yield a unique quantity that is comparable for different values of K, that we denote BIC henceforth.

The BIC is theoretically asymptotically optimal for model selection purpose [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF], however, it may fail in practice for several reasons. In particular, it relies on some hypotheses for the data, such as the i.i.d structure of the residuals, which is violated in fMRI. This means that the goodness of fit of over-parametrized models increases faster than it should in theory, and thus that more complex models, i.e. with a large number of parcels, are systematically and spuriously preferred. In the case of brain volume parcellation, the violation of the i.i.d. hypothesis might be related to different factors, such as data smoothness or spatial jitter across individuals.

Cross-validated likelihood A nice feature of the model ( 7) is that it can be evaluated on test data, thus making it possible to run a cross-validation procedure on different subjects; such a procedure does not overfit, where overfit means models non-reproducible noise, creating the optimistic bias inherent when learning and evaluating a model on the same data. We use the log-likelihood in a shuffle-split crossvalidation scheme: for each fold, the model is learned on the training set (i.e. a random subsample of 80 % of the data): this includes the estimation of the clustering and fitting the mixed-effects model; the log-likelihood computed on the test data is then summed across parcels in order to yield a unique quantity, denoted CV -LL in the following.

Reproducibility by bootstrap

The two previous metrics only address the fit of the data by the model.

Another important criterion in neuroimaging is reproducibility [START_REF] Laconte | The evaluation of preprocessing choices in single-subject BOLD fMRI using NPAIRS performance metrics[END_REF], which we define in this context as the consistency of two clustering solutions across repeats on bootstrap samples taken from the data, measured by assignment statistics of voxels to clusters. To estimate reproducibility, we repeated the clustering by bootstrapping over subjects and assessed the stability of the clustering between pairs of bootstrap samples using two standard metrics: adjusted mutual information or adjusted rand index.

The adjusted Rand index (ARI) is comprised between -1 and 1, and measures the consistency of the two labellings while being invariant to a permutation of the labels [START_REF] Vinh | Information theoretic measures for clusterings comparison: is a correction for chance necessary?[END_REF]. A value of 1 means perfect correspondence of the labeling, while a value of 0 implies that the correspondence are at chance.

An important feature of the ARI metric is that it scales well when the number of clusters K is large. See http://en.wikipedia.org/wiki/Rand_index for more details.

Adjusted mutual information (AMI) upper bounded by 1, and possibly negative, is an estimate of the mutual information of two discrete assignment of voxels to parcels, which is corrected for chance: two statistically independent assignments should have an AMI value of 0, while two identical assignments should have an AMI value of 1 http://en.wikipedia.org/wiki/Adjusted_ mutual_information.

Template subject 0 subject 1 subject 2 subject 3 subject 4 subject 5 subject 6 subject 7 subject 8 subject 9

Figure 1. Example of simulated data used in the 2D simulation experiment. The template or ground truth labeling is shown on the left side, and 10 individual datasets are sampled according to the model, jittered spatially by 2 pixels and then smoothed with a kernel of fwhm 1.17 pixels.

IMPLEMENTATION

We use the algorithms and metrics from the scikit-learn toolbox [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. In particular,

EXPERIMENTS

SIMULATED DATA

Data are simulated according to model ( 6): on a 2D grid of shape 20× 25 pixels, 5 random clusters are generated with a hierarchical clustering approach, by using Ward's parcellation on a set of random signals; 10 individual datasets are sampled using the generative model: for each parcel the parameters µ are sampled from N (0, 1), σ 1 = 1 and the random subject effect β are drawn from N (0, 1), σ 2 = 1. Note that the βs are kept constant across parcels. Data corresponding to a sample of 10 subjects are generated.

To make the data more realistic, we add a deformation to each individual dataset that has a magnitude of 0, 1 or 2 pixels in each direction and smooth it -or not-with a kernel of full width at half maximum (fwhm) of 1.17 pixel. Note however that this breaks (on purpose) the hypotheses of the generative model and makes the simulations more realistic. An example is shown in Fig. 1.

The question that we address is whether we can hope to recover the true number of clusters from the simulation; to do so, we can use one of the three selection criteria: BIC, cross-validation and bootstrap reproducibility (we use B-AMI by default, but B-ARI yields similar results on this dataset). The recovery is quantified through the adjusted rand index between the true labeling of voxels and the obtained one. The results are based on 200 replications of the experiment, and the optimal number K of parcels is searched in the {2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 30} set.

FUNCTIONAL LOCALIZER DATA

Data were acquired from 128 subjects who performed a functional localizer protocol as described in Pinel The data was subject to a pre-processing procedure that includes the correction of the difference in slice timing, motion estimation and correction, co-registration of the EPI volumes to the T1 image, non-linear spatial normalization of T1 images, then of the fMRI scans to the SPM T1 template. All of these steps were performed using the SPM8 software. Optionally, we considered a 5mm isotropic smoothing the normalized images. In parallel, an average mask of the gray matter was obtained from the individual normalized anatomies, subsampled at the fMRI resolution, and used to mask the volume of interest in the functional dataset. This procedure yields approximately Q = 57, 000 voxels at 3mm resolution.

A General Linear Model (GLM) analysis was applied for the volume using the Nipy package http: //nipy.sourceforge.net/. The model included the ten conditions of the experiments convolved with a standard hemodynamic filter and its time derivative, a high-pass filter (cutoff:128s); the procedure included an estimation of the noise auto-correlation using an AR(1) model.

Activation maps were derived for six functional contrasts, that display the activations related to left versus right button presses, motor versus non-motor tasks, sentence listening versus sentence reading, computation versus sentence reading, reading versus passive checkerboard viewing, vertical versus horizontal checkerboard viewing. We consider that these six contrasts give the most usable summary of the topographic information conveyed by the initial ten conditions, without obvious redundancies, while avoiding non-specific effects.

The standardized effects related to these F = 6 contrasts are used for parcellation fit and evaluation. We consider the possible range of values for K: 10, 20, 30,40,50,70,100,150,200,300,400,500,700,1000,1500,2000,3000,5000,7000,10000. We consider the value of the different criteria for different K values.

HCP DATA

A set of N = 67 subjects of the Human Connectome Project (HCP) dataset was also used in our experiments. These subjects are part of the Q2 release; we used the task-fMRI dataset, that comprises 7 different sessions (see [START_REF] Barch | Function in the human connectome: task-fmri and individual differences in behavior[END_REF] for details), all of which are used here. Starting from the preprocessed volume data provided by the HCP consortium, these dataset were analyzed similarly to the Localizer dataset, using the Nipy software for the GLM analysis, that was carried out using the paradigm information provided with the data. The same gray matter mask was used as for the Localizer dataset was used to facilitate comparisons between the two datasets.

In order to reduce computation time, a subset of F = 9 functional contrasts were used: the faces-shape contrast of the emotional protocol, the punish-reward contrast of the gambling protocol, the math-story contrast of the language protocol, the left foot-average and left hand-average contrasts of the motor protocol, the match-relation contrast of the relational protocol, the theory of mind-random contrast of the social protocol and the two back-zero back contrast of the working memory protocol. this choice was meant to sample a significant set of cognitive dimensions tested in the protocol, without being exhaustive.

This is a provisional file, not the final typeset article [2,3,4,5,6,7,8,9,10,15,20,30] and that the true value is 5. The results are based on data smoothed with a kernel of size 0.5 voxel, and under spatial jitter of 1 voxel isotropic, and are presented across 200 replications. bic, cv (cross validation) and repro (Adjusted Mutual Information) represent three model selection approaches, while ward, k-means and spectral represent three different clustering approaches. After parcellation, the parameters of the model 6 are estimated in each parcel, for each functional contrast and can be plotted in the brain volume; see Fig. 4. In particular, it can be seen that σ 1 > σ 2 uniformly i.e. within-parcel variability dominates across-subject variability when K = 500. Moreover, in the case of Ward's parcellation presented here, the within-and between-subject variance estimated are quite homogeneous across the brain volume. Note however that there is a tendency for both to be correlated with the absolute value of the mean signal. Next, we consider how the variance components, averaged across parcels, change with K in Fig. 5. These values evolve monotonously with K: the intrasubject parameter σ 1 (that measures the cross-voxel variance within a given subject, averaged across parcels) decreases monotonously with K, as expected; the inter-subject parameter σ 2 , that characterizes the cross-subject variability of the mean signal within a parcel, increases monotonously. Both parameters come close to equality for large values of K (about 5000). These trends are similar across clustering techniques. This actually means that changing the resolution yields a re-allocation of the variance from the intra-subject to the inter-subject component of the mixed-effects model. More specifically, for low values of K, the high within-subject variance shadows the between-subject variance, and a very large value of K has to be used if one wants to estimate correctly the between-subject variability of the BOLD signal within the parcellation framework. . Dependence on K of the variance components from model 6, averaged across parcels and contrasts: both σ 1 and σ 2 parameters show a monotonic behavior: the within subject variance decreases σ 1 with K, while the between-subject variance σ 2 increases with K.

RESULTS

SIMULATIONS

Atlas

Summed log-likelihood stdv Harvard-Oxford atlas -6.642 10 7 1.9 10 5 Geometric parcellation -6.589 10 7 1.9 10 5 k-means parcellation -6.463 10 7 1.8 10 5 Ward parcellation -6.513 10 7 1.9 10 5 Spectral clustering parcellation -6.591 10 7 1.8 10 5 (HCP) k-means parcellation -6.710 10 7 1.9 10 5 (HCP) Ward parcellation -6.522 10 7 1.9 10 5 (HCP) Spectral clustering parcellation -6.613 10 7 1.9 10 5

Table 1. Summed log-likelihood of the Localizer data under different spatial models (the higher, the better): brain atlas (top), parcellation on the Localizer dataset (middle), parcellations from the HCP data (right). The number of parcels used is K = 158 for all methods. The standard deviation is obtained by drawing 30 bootstrap samples.

Comparison with an anatomical atlas As a basis for comparison with anatomical atlases, we evaluated the log-likelihood of the data with the most detailed atlas that we could find. We used the Harvard-Oxford atlas both cortical and subcortical http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Atlases together with the cerebellar atlas [START_REF] Diedrichsen | A probabilistic mr atlas of the human cerebellum[END_REF]. The version used was that of FSL 4.1. The regions were systematically divided into left and right hemispheres by taking the sign of the x MNI coordinate of the voxels. Using this procedure, we obtained 158 regions. This atlas was resampled at the resolution of the test fMRI data, and the likelihood of the data summed over parcels was evaluated and compared with that of data-driven parcellations with 158 parcels, obtained either from the Localizer dataset itself or from the HCP dataset. Standard deviation of the log-likelihood are obtained by drawing B = 30 bootstrap samples. The results are shown in Table 1.

We show the corresponding results on the HCP dataset in Table 2.

It can be seen that the anatomical atlas achieves the poorest fit: summarizing the fMRI data on the corresponding set of parcels looses a lot of information. Even a purely geometric parcellation performs better, which can be understood given that it tends to create parcels with equal size, hence achieves a more Table 2. Summed log-likelihood of the HCP data under different spatial models (the higher, the better): brain atlas (top), parcellation on the HCP dataset (middle), parcellations from the Localizer data (right). The number of parcels used is K = 158 for all methods. The standard deviation is obtained by drawing 30 bootstrap samples.

regular sampling of the volume of interest. For K = 158 the best performing parcellation on the training set is obtained from k-means, but these parcellations do not generalize well from a dataset to another.

Ward's parcellation on the other hand, performs better than geometric clustering in all configurations.

Finally, the bootstrap variability of these results is typically small with respect to between-method difference for the Localizer dataset, ensuring that the differences are significant. This is less so with the HCP dataset, for two reasons: the number of subjects is smaller, and the per-subject SNR seems relatively lower in that dataset (see [START_REF] Barch | Function in the human connectome: task-fmri and individual differences in behavior[END_REF]).

Analysis of the goodness of fit the models (Localizer dataset)

The goodness of fit of the model is given by the log-likelihood, which can be compared across methods for a fixed value of K in Fig. 6 (a).

The main observations are:

• For all methods, the curve achieves an optimum value for a very large number of parcels (3000 ≤ K ≤ 7000), which is much more that the number typically expected and used in neuroimaging experiments.

• k-means and Ward's clustering achieve the lowest distortion -i.e. loss of information from the original signal-, with k-means performing better for low number of parcels and Ward's clustering performing better for large number of clusters. Spectral clustering is inferior in terms of goodness of fit. It is even lower than a purely geometric parcellation of the brain volume for some values of K.

• The achieved log-likelihood is larger on smoothed data than on unsmoothed data, but the behavior is qualitatively similar. In this report, we present only results on unsmoothed data.

Second, we can observe that, unlike in our simulations, BIC and cross-validated log-likelihood (Fig. 6 (b-c)) achieve their optimum at the same value of K as the data log-likelihood function, thus at very high values (3000 ≤ K ≤ 7000). This late increase can readily be interpreted as an artifact due to the fact that we are now observing a very 6.9 6.8 6.7 6.6 6.5 6.4 6.3 6.2 6.1

Log-likelihood large number of very small clusters, and that the reproducibility indexes are not well suited in this case.

It is also true that very small clusters tend to represent the spatial neighboring system, and thus this high reproducibility is not very informative on the functional features carried by the data.

The spectral clustering outperforms the other alternatives regarding reproducibility, which means that it is able to capture some stable features in the input data, although the overall representation is suboptimal in terms of accuracy. Regarding the sensitivity/reproducibility compromise (see Fig. 7, right), the spectral method is dominant in the low accuracy/high reproducibility region, while Ward's method dominates in the high accuracy/low reproducibility region.

Frontiers in Brain Imaging Methods

MODEL SELECTION RESULTS ON THE HCP DATASET

A summary of the results obtained by doing the same experiments on the HCP dataset is provided in Fig. 8. In spite of weak changes of the optimal values K ⋆ ≤ 3000 for accuracy, K ⋆ ∈ [200, 500] for reproducibility), this dataset reproduces exactly the trends observed with the Localizer dataset: Ward's method outperforms the others in terms of accuracy and for high K values, spectral clustering yields a poor fit and a high reproducibility, k-means a good fit, especially for small K, yet very low reproducibility.

DISCUSSION

Our experiments benchmark three methods to derive brain parcellations from functional data, using three model selection criteria and two reproducibility measures. Though not exhaustive, these experiments are very informative on the general behavior, the domain of optimality of the methods, and the issues that limit the power of such approaches in neuroimaging data analysis.

GUIDELINES FOR FUNCTIONAL PARCELLATION EXTRACTION

Which criterion to use for methods comparison? To frame the problem, it is necessary to choose the criterion used to guide model selection. Note that this is an important yet difficult aspect of any unsupervised statistical learning procedure. We studied two different characteristics of functional parcellation that are critical to their usage in brain mapping: how well they capture the functional signal and how reproducible they are under perturbations of the data. To measure the goodness of fit of the functional signal, it is important to distinguish within-subject variance from across-subject variance, as only the first kind of variance is minimized when the number of parcels increases. Our probabilistic model offers a natural goodness of fit criterion, the log-likelihood; by penalizing it (BIC criterion) or by using cross-validation, it is possible to obtain a sound model selection. Our simulations show that cross-validation almost systematically outperforms BIC, but we did not notice systematic differences in the real dataset. The other important aspect of a brain description is its stability, and we also investigated other criteria that the selected number of clusters based on the consistency of parcellations. This approach behaved similarly as the others on synthetic data, but provided a much more conservative selection on real data (K ∼ 200 to 500 parcels according to the dataset and method). The fact that reproducibility and accuracy yield different decisions for model selection is well known, and has been illustrated in neuroimaging by [START_REF] Laconte | The evaluation of preprocessing choices in single-subject BOLD fMRI using NPAIRS performance metrics[END_REF]; this effect is tightly related to the classical bias /variance in statistics.

Which algorithm to prefer? Regarding the clustering algorithms themselves, our general finding is that Ward's algorithm should be preferred, unless a small number of parcels is required. Indeed, spatiallyconstrained Ward's clustering outperforms the other techniques in the large K regime (say, K ≥ 500)

in terms of goodness of fit, while having fair results in terms of reproducibility. With respect to kmeans, it offers the additional advantage of providing spatially connected parcels. In theory, k-means algorithm should do better in terms of accuracy, but the optimization problem solved by k-means is hard (non-convex) and thus bound to sub-optimal solutions; as a consequence the greedy approach in Ward's algorithm outperforms it. Moreover, k-means based parcellations tend to fit data idiosyncrasies and thus do not generalize well across datasets, as shown in tables 1 and 2. We observed that mixture models would behave similarly to k-means, since k-means is in fact a constrained Gaussian mixture model with hard assignments. In a side experiment, we observed that Gaussian mixture models perform consistently better than k-means, but the difference is tiny and comes with high computational cost.

Spectral clustering is not a powerful approach to outline structures in the data. A simple geometric clustering procedure is as good, and sometimes better in terms of accuracy. The reason is that spectral clustering is efficient with high SNR data when clusters are easily discriminated, which is not the case with functional neuroimaging data, where it mostly outlines geometrical structures. Similar observations were made in [START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF]. Note however that spectral clustering is even more stable than geometric clustering, meaning that it captures some structure of the input data.

How many parcels? It should first be emphasized that choosing the number of parcels in our model is not exactly the question of deciding how many functional regions can be found in the brain, but how many piecewise constant models can actually be fit to some fMRI data reliably. The distinction is important, because some regions, for instance V1, will contain internal functional gradients, such as those related to retinotopy, orientation sensitivity and ocular dominance. In theory, the function specificity could therefore be resolved at the level of columns in these regions, but this does not mean that larger structures do not exist. The conclusions that we draw here are bound to the data that we have used and generalization to different modalities or contrasts (resting-state fMRI, anatomical connectivity) is not guaranteed.

The goodness of fit-related criteria yields high numbers (up to K = 5000 for Ward's clustering, slightly less for the others, but this may simply reflect a lack of sensitivity of these approaches in the large K regime, in which Ward's clustering fits the data better), simply indicating that functional activations cannot easily be represented as piecewise constant models: whether this is an intrinsic feature of brain function or an impact of cross-subject spatial mismatch or pre-processing artifacts remains an open question. In the future, the use of brain registration algorithms based on functional data [START_REF] Sabuncu | Function-based intersubject alignment of human cortical anatomy[END_REF]

  et al. (2007) and referred to as Localizer henceforth. This protocol is intended to activate multiple brain regions in a relatively short time (128 brain volumes acquired in 5 minutes) with ten experimental conditions, allowing the computation of many different functional contrasts: left and right button presses after auditory or visual instruction, mental computation after auditory or visual instruction, sentence listening or reading, passive viewing of horizontal and vertical checkerboards. The subjects gave informed consent and the protocol was approved by the local ethics committee.In 59 of the subjects, functional images were acquired on an 3T Siemens Trio scanner using an EPI sequence (TR = 2400ms, TE = 60ms, matrix size = 64 × 64, FOV = 19.2cm × 19.2cm). Each volume consisted of 40 3mm-thick axial slices without gap. A session comprised 132 EPI scans, of which the first four were discarded to allow the MR signal to reach steady state. The slices were acquired in interleaved ascending order. Anatomical fSPGR T1-weighted images were acquired on the same scanner, with a slice thickness of 1.0 mm, a field of view of 24 cm and an acquisition matrix of 256 × 256 × 128 voxels, resulting in 124 contiguous double-echo slices with voxel dimensions of (1.0 × 1.0 × 1.0) mm 3 .In 69 of the subjects, functional and anatomical were acquired on a 3T Bruker scanner. Functional images were acquired using an EPI sequence (TR = 2400ms, TE = 60ms, matrix size[START_REF] Golland | Detection of spatial activation patterns as unsupervised segmentation of fmri data[END_REF] FOV = 19.2cm × 19.2cm). Each volume consisted of n a 3-mm-or 4-mm-thick axial slices without gap, where n a varied from 26 to 40 according to the session. A session comprised 130 scans. The first four functional scans were discarded to allow the MR signal to reach steady state. Anatomical T1 images were acquired on the same scanner, with a spatial resolution of (1. × 1. × 1.2) mm 3 .

Figure 2 .

 2 Figure2. Results of the simulations: choice of the number parcels for different clustering methods and cluster selection techniques. Note that the range of possible values is[2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 30] and that the true value is 5. The results are based on data smoothed with a kernel of size 0.5 voxel, and under spatial jitter of 1 voxel isotropic, and are presented across 200 replications. bic, cv (cross validation) and repro (Adjusted Mutual Information) represent three model selection approaches, while ward, k-means and spectral represent three different clustering approaches.

Fig. 2

 2 Fig. 2 displays the selected K ⋆ value, based on data smoothed with a kernel of size 0.5 voxel, and under 298

Figure 3 .

 3 Figure 3. Example of parcellation with 500 parcels on the Localizer dataset.

Figure 4 .

 4 Figure 4. Example of parameters estimated in a parcellation obtained with Ward's clustering and K = 500 parcels. They are given in arbitrary units (percent of the baseline fMRI signal, squared for variance estimates). These parameters are those for the computation-sentence reading functional contrast.

Figure 5

 5 Figure5. Dependence on K of the variance components from model 6, averaged across parcels and contrasts: both σ 1 and σ 2 parameters show a monotonic behavior: the within subject variance decreases σ 1 with K, while the between-subject variance σ 2 increases with K.

  -reproducibility compromise (Localizer dataset) The reproducibility of the clustering estimated by bootstrapping the data can be studied as a function of the number of clusters, or as a function of the likelihood. Both representations are presented in Fig. 7.The reproducibility index displays a clear optimum value at K ≃ 200 parcels. For larger values, the reproducibility index decreases slowly, but increases again for very large number of parcels K > 4000.

Figure 6 .

 6 Figure 6. Distortion metrics (a) Accuracy of the model 6 measured through the summed Log-likelihood across parcels, as a function of the number K of clusters. The accuracy is maximized for very high values of K. The Bayesian Information Criterion (b) -with the sign flipped for the sake of visualization-and the cross-validated log-likelihood (c), that can be used to identify the right model show the same behavior as the log-likelihood function.

Figure 7 .

 7 Figure7. Analysis of the reproducibility index with respect to the number of parcels (a, b) and with respect to the negative log-likelihood (c). For all methods but one, the B-AMI (Bootstrapped Adjusted Mutual information) index (a) shows a (local) maximum for about 200 parcels and decreases against for larger numbers, until it increases again for very large number clusters (K ≥ 5000). By contrast, B-ARI (Bootstrapped Adjusted Rand Index) (b) only displays the local maximum on Ward's parcellation. If we consider B-AMI against accuracy there is thus a trade-off region, for a number of parcels comprised between 200 and 5000 (decreasing portion of the curves in the reproducibility-accuracy curve), in which each setting represents a different compromise. The two dominant techniques are spectral clustering, that maximizes the reproducibility index, and Ward's clustering, that yields higher accuracy overall.

Figure 8 .

 8 Figure 8. Results of the of the model selection experiments on the HCP dataset: (left) accuracy-based selection through the BIC score, (middle) reproducibility-based selection through Bootstrapped Adjusted Mutual Information, (right) ensuing sensitivity/reproducibility curve.

  may significantly affect model selection. The reproducibility criterion, on the other hand, peaked at K ∼ 200, meaning that there is probably a relevant level of description with such a resolution. Thus, when parcellations are used to obtain a model of brain organization that seeks to characterize individually each parcel, a conservative choice K ∼ 200 to 500 should be preferred for the sake of reproducibility. Note that K = 200 is a lower bound on the right dimensionality, i.e. models with a resolution lower than 200 regions are not flexible enough to represent functional signals without introducing severe distortions. In particular, anatomical atlases that propose a decomposition into about 100 regions, are not sufficient to summarize functional signals, some of the resulting ROIs being very large. Yet, the problem of optimizing the number of parcels remains open and should be addressed in a datadriven fashion.

  Ward's algorithm is very efficient on data size of typical brain images. The following version of the software were used:Matlab R2013A, version 8.1.0.64, SPM8 v. 5242, scikit learn v. 0.14. The code used in this work is available at https://github.com/bthirion/frontiers_2014.

Frontiers in Brain Imaging Methods

This is a provisional file, not the final typeset article

A variant of spectral clustering replaces this k-means step by learning a rotation to discretize the representation[START_REF] Yu | Multiclass spectral clustering[END_REF]. We used this approach, that outperforms k-means, in the case where the number of desired clusters K is smaller than the subspace dimension m (see next).Frontiers in Brain Imaging Methods

ACKNOWLEDGMENT

We would like to thank Philippe Pinel for providing the Localizer data. The so-called HCP data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University Funding: The authors acknowledge support from the ANR grant BrainPedia ANR-2010-JCJC-1408-01 and the Human Brain Project.

CHALLENGES AND FURTHER WORK

The difficulty of model selection on noisy data It is important to remember that discovering functionally homogeneous structures is a hard problem, given that the SNR of the data is low, and that even visual inspection would most often be insufficient to define a relevant parcellation. Besides this issue, neuroimaging data come with additional difficulties: the data are smooth, which could be accounted for but is not in model 6. The other difficulty is that the spatial jitter brought by imperfect spatial normalization, and poor matching between functional organization and sulco-gyral anatomy across subjects makes this an ill-posed problem, since regions with homogeneous functional characteristics may be slightly displaced across individuals, which invalidates the model hypotheses. Both smoothing and jitter break the hypotheses of BIC, which yields poor model selection. Cross-validation and reproducibility are more resilient to this effect.

Limitations of this experiment Our experiments are based on two datasets, with a pre-defined set of contrasts. We have been able to check that using any of the contrasts or all of them yields qualitatively similar results (data not shown). The power of the experiment is that it is based on a relatively large number of subjects (67 to 128), so that one can at least conjecture that the between-subject variability observed in functional neuroimaging is correctly sampled. Note that the Localizer data come from two different scanners, resulting in an un-modeled latent factor. We observed, however, that our conclusions were unaltered when performed on a subset of subjects coming from the same scanner (data not shown).

The model that we use has several limitations:

• The parcellation itself is fixed across subjects. While a relaxation to individual dataset has been proposed in [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: multi-subject parcellation of fmri datasets[END_REF], such a procedure loses some of the properties of clustering, and make model selection much harder.

• Our model (eq. 6) does not account for spatial effects in the within-parcel covariance, which would probably make it more robust to data smoothness and possibly to cross-subject spatial jitter, but the computational price to pay for these models is high.

• It assumes that the true activation signal is piecewise constant. A smooth interpolation scheme between parcels might make it more powerful, hence reducing the requirement of large K values.

Again, this would increase the complexity of the model fitting.

Suggestions for population-level fMRI modeling One of the observations made in this study is that the problem of the spatial jitter across subjects remains the main limitation that needs to be overcome in order to learn appropriate population-level atlases. This should be addressed using procedures such as those presented in [START_REF] Sabuncu | Function-based intersubject alignment of human cortical anatomy[END_REF]; [START_REF] Robinson | Multimodal surface matching: Fast and generalisable cortical registration using discrete optimisation[END_REF]. Other improvements of the model concern the possibility of using not a single parcellation, but several different parcellations and to aggregate the results (i.e. the significant effects across subjects) by marginalizing the parcellation as a hidden variable of parametric models [START_REF] Da Mota | Enhancing the Reproducibility of Group Analysis with Randomized Brain Parcellations[END_REF][START_REF] Varoquaux | Small-sample brain mapping: sparse recovery on spatially correlated designs with randomization and clustering[END_REF]. Besides, different parcellation schemes could use different values for K. In particular, Ward's algorithm is a hierarchical algorithm, that can actually be used to estimate multi-scale representations of brain activity (see e.g. [START_REF] Michel | A supervised clustering approach for fMRI-based inference of brain states[END_REF]; [START_REF] Orban | The richness of taskevoked hemodynamic responses defines a pseudohierarchy of functionally meaningful brain networks[END_REF]). Specifically [START_REF] Orban | The richness of taskevoked hemodynamic responses defines a pseudohierarchy of functionally meaningful brain networks[END_REF] suggest that the hierarchical organization of nested clusterings obtained from hemodyamic response function would be stable in the population, hinting at an intrinsic feature of brain organization. This is an additional asset of this procedure that has not been considered in this work but could be used in future applications of brain parcellations.

DISCLOSURE/CONFLICT-OF-INTEREST STATEMENT

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


Bertrand Thirion

Gaël Varoquaux

Olivier Grisel

Cyril Poupon

Philippe Pinel

Principal Component Regression predicts functional responses across individuals

Keywords: fMRI, principal components regression, random effects

Inter-subject variability is a major hurdle for neuroimaging group-level inference, as it creates complex image patterns that are not captured by standard analysis models and jeopardizes the sensitivity of statistical procedures. A solution to this problem is to model random subjects effects by using the redundant information conveyed by multiple imaging contrasts. In this paper, we introduce a novel analysis framework, where we estimate the amount of variance that is fit by a random effects subspace learned on other images; we show that a principal component regression estimator outperforms other regression models and that it fits a significant proportion (10% to 25%) of the betweensubject variability. This proves for the first time that the accumulation of contrasts in each individual can provide the basis for more sensitive neuroimaging group analyzes.

Introduction

Functional brain imaging is a tool of choice to understand the functional specialization of brain areas and decipher the mechanisms or the impact of various brain diseases. It outputs contrast maps that quantify metabolic responses of the brain gray matter in various experimental situations; more precisely, a contrast is a certain combination of the activation evoked by various experimental conditions that reveals brain responses specific to a cognitive process; the accumulation of such contrasts can then be used to derive subject-specific or population-level functional cartography [START_REF] Lashkari | Search for patterns of functional specificity in the brain: a nonparametric hierarchical bayesian model for group fMRI data[END_REF]. The underlying inference is based on a mass-univariate framework that detects which brain voxels show in average a positive response. Between-subjects variability enters the model as a random effect: larger variability entails lower evidence that an effect is present, hence less sensitivity [START_REF] Worsley | A general statistical analysis for fMRI data[END_REF].

While classical analyzes ignore the structure underlying these random effects, there is a growing evidence that some of it can be characterized or even predicted by several factors, such as i) the presence of an average effect in the population (high variability coincides spatially with larger mean effect) [START_REF] Thirion | Analysis of a large fMRI cohort: Statistical and methodological issues for group analyses[END_REF] ii) anatomical differences, such as the cortical ends of anatomical connections, which has been illustrated by a prediction of faces versus scenes contrast in the fusiform gyrus [START_REF] Saygin | Anatomical connectivity patterns predict face selectivity in the fusiform gyrus[END_REF] and iii) individual functional connectivity between brain regions, which has been shown to provide a relevant prior for group analyzes in [START_REF] Ng | Connectivityinformed fMRI activation detection[END_REF]. Modeling such random effects is important in order to facilitate future correlative studies that compare brain features to genetics and behavior.

In this work, we consider the prediction of functional contrasts by other contrasts observed independently in the same subjects, assuming that the betweensubject variability can be captured through a latent subspace structure. This approach is reminiscent of the popular modality propagation analysis paradigm (among others, see [START_REF] Ye | Modality propagation: coherent synthesis of subject-specific scans with data-driven regularization[END_REF]), in which a target image is reconstructed from a reference one by machine learning tools; the main difference is that functional neuroimaging data have a much lower signal-to-noise ratio. However, a large number of contrast images can potentially be used to capture the structure of the random effects space. Altogether, this implies that the same procedures are not optimal for random effects characterization and modality prediction.

The contribution of this paper is two-fold: i) we propose a formalization of the random effects subspace framework and derive effective methods to estimate such models on real data; in particular, principal component regression (PCR [START_REF] Dhillon | A risk comparison of ordinary least squares vs ridge regression[END_REF]) comes naturally as a tool of choice in this problem; ii) we benchmark several estimators to test the effectiveness of PCR in that setting as well as three spatial frameworks in the inference process: a local framework based on each spatial unit of the domain, a regional model that uses information from neighboring regions, and a global framework that estimates the random effects subspace on an initial decomposition of the images into principal components. [START_REF] Worsley | A general statistical analysis for fMRI data[END_REF] The random effects subspace model Let (X s ), s ∈ {1, .., n subjects } be (n contrasts , n voxels ) matrices that represent a set of imaging contrasts observed in a given group of subjects. X s is typically analyzed as the sum of a fixed effect X across individuals, random effects that represent the between subject variability and observation noise. The simplest approach to extract the latent structure of these random effects is a factorization that represents X s -X as a combination of subject-specific latent factors ζ s = (ζ k ) s , k = 1..K of shape (n voxels ), weighted by voxel-specific coefficients α = (α k ), k = 1..K that are common to to the population:

X s = X + K k=1 ζ s k α k + ε, ∀s ∈ {1, .., n subjects } ( 1
)
where ε is the observation noise, K is the number of latent components, and the product ζ s k α k is taken element-wise. Note that both α and ζ s can vary spatially. The identification of model (1) yields a series of questions: i) what fraction of the variance is fit by the latent components versus observation noise? ii) to which extent is the amount of explained variance correlated to the fixed effects? iii) Are the latent components shared across brain regions or specific to a given location? These questions can be handled by predictive models that quantify the amount of residual variance var(X t |X ref ), where X t , X ref correspond respectively to a test and a set of reference contrasts acquired independently in the same subjects. The corresponding model is the following one:

X t s = Xt + K k=1 (ζ ref ) s k α tk + ε, ∀s ∈ {1, .., n subjects } (2)
and it goodness of fit is can be measured by

n subjects s=1 X t s -Xt + K k=1 (ζ ref ) s k α tk 2 .
In this work, we explore several solutions to estimate the latent factors ζ ref :

- -Lasso regression, that performs a non-linear, sparsity-inducing shrinkage on the loadings (α t ), given

ζ s ref = X ref s -Xref .
-Extremely Randomized Trees (Extra Trees) [START_REF] Geurts | Extremely randomized trees[END_REF], that also take

ζ s ref = X ref s -
Xref , yet learn (α t ) by tree-based regression as in random forests.

Next, to address the question whether the latent components are local to some image region, shared locally, or across brain regions, we perform the analysis at three different spatial scales:

-The local scale, in which the model ( 2) is estimated independently at each brain location considered. -The regional scale, in which ζ ref and (α t ) are estimated in the neighborhood of a given image location. -The global scale, in which the model is estimated on a global model of the data. In this work, we use a PCA of the whole dataset, written as a (n voxels , n subjects × n contrasts ) matrix; we used 400 components.

Comparing the predictive scores of these approaches captures a bias/variance trade-off -local estimators being less biased but more variable-but also indicates how redundant the information is, as redundancy is used only by global or regional methods. The learners used in this paper are the scikit-learn implementation [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] v0.14, used with default parameters. The code can be found at https://github.com/bthirion/fMRI_PCR.

Experiments

As our aim is model comparison, we directly turn to an actual MRI dataset to assess how well test contrasts are predicted by the above estimators. Data. A dataset consisting of four functional MRI sessions acquired under different experimental protocols and an anatomical image, obtained in 79 subjects, was used. Standard preprocessing were carried out with SPM8 software and its default parameters; functional images were warped in the MNI152 coordinate space and resampled at 3mm resolution. A population mask of the gray matter voxels was obtained by averaging individual gray matter probability maps obtained from SPM8, resampling this average at the fMRI data resolution and thresholding it, which yields about 57,000 voxels. Contrasts were obtained using a standard linear model, based on the specified task timing information,canonical hemodynamic response function, high-pass filtering and auto-regressive noise model. The estimation of the model parameters was carried out using the Nipy software.

Contrasts. The four functional protocols were the following: i) a generic localizer protocol, that assesses basic functions, such as finger tapping, listening and reading, attentional tasks; ii) a social protocol in which the subjects were exposed to false beliefs, to objects moving with or without a putative intention, speech and non-speech iii) a so-called parietal mapping protocol, during which the subjects had to perform visual saccades, grasping, judge the orientation of an object; iv) an emotional protocol, during which the subjects had to judge the confidence inspired by or the gender of a face or of an expression, together with a control tasks on the same stimuli. Five to six main contrasts were retained from each experiment to form a set of functional contrasts of interest. The reference set consisted of the contrasts i), ii) and iii) or i), ii) and iv), while the test set was the set of contrasts of the left-out acquisition.

Evaluation metric. Based on the approach described in section 2, we use the following criterion:

ρ = 1 -Xt-( Xt|Xref ) 2 Xt-Xt 2
, where the norm is taken over the whole image, computed over a cross-validation loop on subjects. ρ can readily be interpreted as the explained variance ratio of the model (the higher the better). In all experiments presented here, the test contrast has been observed in another session than the reference contrasts, in order to enforce independence of the samples. In our experiments, we used a 20-fold shuffle-split procedure, where 90% of the subjects are used for training and 10% are used for testing. Parameter setting. All the free parameters of the methods (amount of penalty for the ridge and lasso models, number of principal components for PCR) were set by nested cross-validation, using 3-fold cross-validation on the training set. Finally, the data presented here are based upon an initial parcellation of the dataset into Q = 10, 000 super-voxels based on Ward clustering instead of the 57, 000 voxels of the reference cortical mask. Voxel-specific contrasts were averaged to parcelspecific values. Note that this procedure is simply meant to be a compression of the data and has no influence on the model nor its interpretation. We checked that using finer or coarser parcellations does not alter the experimental results.

Results

Method comparison. A boxplot of the ρ score computed over 11 contrast is presented in Fig. 1. Regarding estimator comparison, this shows that, for all spatial scales under consideration (local, regional, global), the following hierarchy holds: PCR outperforms Lasso, that outperforms Extra-Trees regression and the ridge estimator. These differences significant using pairwise Wilcoxon tests (PCR >lasso, lasso > ridge, lasso > extra trees, p < .005, in all spatial configurations; Paired Student tests yield much lower p-values). Note that we considered other variants of ensemble estimators, such as Gradient Boosted Trees and Random Forests, but these were systematically outperformed by the Extra Trees estimator. Second, local estimators significantly (p < 0.005) outperform regional and global estimators: this means that there is not enough information shared across regions to improve the predictive accuracy of the model, while the addition of new regressors degrades its predictive accuracy, through a classical curse of dimensionality phenomenon.

Well-predicted brain regions are those that show large fixed effects. To better understand the mechanisms of across-contrasts random effect prediction, an important question is whether the prediction accuracy is uniform across brain regions, or whether it coincides with regions that show a strong mean effect. In Fig. 3, we show for a contrast displaying the activation evoked by deciphering the intention from face expression, the map of the (across subjects) mean statistical response together with a map of the prediction accuracy obtained with local PCR. This shows that high values are obtained in the same occipital, parietal and pre-frontal cortical regions. This effect holds for all other contrasts: the bottom part of Fig. 3 shows the joint density of ρ versus the average (across contrasts and subjects) contrast absolute value (to discount sign effects); this shows that the regions with high predictability are also those with larger mean activation strength. The 10% of the regions showing the largest average responses and those showing highest predictability (across contrasts and subjects) are displayed in the same figure (right side): these are the occipital, motor cortex and supplementary motor area in both cases; note however that the orbito-frontal cortex shows high activation in average and relatively lower predictability, while the posterior superior temporal gyrus is highly predictable although the mean activation is weaker there. ROI-based analysis. To better illustrate potential differences between predictability and average activation, we display the average functional signal for two nearby ROIs of the supra-marginal gyrus (x, y, z = 65, -27, 23 and 65, -37, 23 mm respectively): both show the same amount of activation across eleven contrasts, but the most anterior one is much more predictable than those of the posterior one. A straightforward explanation is that the anterior region is much more consistently defined in the MNI space as the border of the post-central gyrus, while the posterior one is a region that does not get a consistent definition, and is more variable anatomically and functionally [START_REF] Plaze | where do auditory hallucinations come from?"a brain morphometry study of schizophrenia patients with inner or outer space hallucinations[END_REF][START_REF] Essen | Parcellations and hemispheric asymmetries of human cerebral cortex analyzed on surface-based atlases[END_REF].

Discussion

This paper presents the first successful attempt to predict individual functional responses based on previously observed data from the same subjects, evaluated across all brain regions. Here the prediction is based on the observation of the same contrasts in other subjects and mimics the realistic incomplete dataset situations where some general contrasts have been observed in some subjects while others are missing. The accuracy of the prediction is relatively high: 10% to 25 % of the variance can be fit using the local PCR model. This confirms a posteriori the existence of a latent subspace that spans the variability of individual responses. Moreover, the best performing technique, PCR, is precisely the one that builds on the low-dimensional latent space model of the variability. The potential superiority of PCR against ridge has already been described [START_REF] Dhillon | A risk comparison of ordinary least squares vs ridge regression[END_REF] and in the present case, it confirms that dimension reduction through PCA of a set of reference contrasts is better at separating random effects from noise than the uniform shrinkage of ridge estimators. Ridge estimators completely fail in the situation where more regressors are used, e.g. in the so-called regional model. Interestingly, lasso estimators have an intermediate behavior, showing that the ℓ 1 shrinkage provides a relevant proxy for subspace estimation, although it performs less well than mere dimension reduction. Unlike other settings, tree-based methods do not perform well here, probably because the sample/dimensionality ratio is not high enough. Local models outperform regional or more global models; however, we conjecture that randomized or multiscale models that average the predictions across spatial models could perform better. The regions that show high predictability are mostly those that show high average activation across subjects and they tend to coincide with the visualsensorimotor-auditory network described in [START_REF] Mesmoudi | Resting state networks' corticotopy: the dual intertwined rings architecture[END_REF], with the addition of the attentional network. Most importantly, the fact that predictability is correlated with average activation means that traditional mass-univariate random effects models used in neuroimaging studies are suboptimal, as they handle the random effects as an error term, while they actually represent more complex structures (mismatch in the coregistration, physiological variability) and may be indicative of relevant anatomical and functional differences. A closer analysis reveals some dissociations: some regions show more variability for a given amount of average activation. This is likely because individual anatomical variability has been captured by the model, an hypothesis suggested in [START_REF] Essen | Parcellations and hemispheric asymmetries of human cerebral cortex analyzed on surface-based atlases[END_REF] that had not been confirmed on a functional basis so far. Individual response prediction can thus complement the traditional view of brain mapping and will pave the way toward more sensitive functional neuroimaging group analyzes.

Fig. 1 .

 1 Fig. 1. Comparison of the methods for prediction of test contrasts. The boxplot represents the average rate of explained variance through different methods for eleven different functional contrasts.

u s t in e s s e x p r . in t e n t io n e x p r e s s io n g e n d e r e x p r e s s io n c o n t r o l o b je c t g r a s p in g o b je c t o r ie n t a t iou s t in e s s e x p r . in t e n t io n e x p r e s s io n g e n d e r e x p r e s s io n c o n t r o l o b je c t g r a s p in g o b je c t o r ie n t a t ioFig. 2 .

 2 Fig.2. Predictability reveals other features than average statistical sensitivity. Two close regions of interest have been selected, in the right supra-marginal gyrus (x=65mm, y=-27 or -37mm, z=23mm ). While both regions approximately display similar activation strength across 11 test conditions (left), the individual responses in the anterior one is much more predictable than the responses from the posterior region (right).

Fig. 3 .

 3 Fig. 3. Coincidence of regions that show high activation values and those that are actually predictable across individuals. (Top) Cross-subject average activation level for a contrast related to deciphering expression intention. (Middle) Rate of prediction across the same regions for the same contrast. (Bottom) left: joint histogram of the mean statistical significance and ratio of explained variance per voxel, across experimental conditions; middle: 3D view of the most predictable regions across conditions; right: 3D view of the regions that display the strongest average responses across conditions.

  Principal component regression (PCR), that defines ζ ref by a PCA of X ref and estimates the loadings (α t ) by simple linear regression.

-Ridge regression considers ζ s ref = X ref s -Xref , but applies a uniform shrinkage on the loadings (α t ).

Acknwoledgements The authors acknowledge support from the ANR grant Brain-Pedia ANR-2010-JCJC-1408-01 and the Human Brain Project. They are very thankful to D. lebihan, the PI of the archi study, for sharing the data.