Bertrand Thirion

Philippe Pinel

Alan Tucholka

Alexis Roche

Philippe Ciuciu

Jean-François Mangin

Jean-Baptiste Poline

Structural analysis of fMRI data revisited

Keywords: functional MRI, Group analysis, spatial normalization, structural methods, watershed, belief propagation, replicator dynamics, group comparison

Structural Analysis of fMRI Data Revisited:

Improving the Sensitivity and Reliability of fMRI Group Studies I. INTRODUCTION Functional neuroimaging aims at finding brain regions specifically involved in the performance of cognitive tasks. In particular, functional MRI (fMRI) is based on the detection of task-related Blood Oxygen-Level Dependent (BOLD) effect in the brain. The measurement of this effect is performed by regression analysis of four-dimensional datasets (three spatial dimensions plus time) against pre-defined regressors that represent the expected BOLD response to the stimulations across time; this analysis framework is known as the General Linear Model (GLM) [START_REF]Human Brain Function[END_REF]. Inference about putative regions of activity is generally based on several subjects (∼10-15 subjects typically), and the current standard procedure consists in detecting voxels for which the average task-related BOLD signal increase is significant across subjects (random/mixed effects analyses, R/MFX) [START_REF]Human Brain Function[END_REF], [START_REF] Worsley | A general statistical analysis for fMRI data[END_REF]. Such voxel-based inference schemes require the images to be warped to a common space, which is usually performed by coregistration of the anatomical, then functional data with a template image [START_REF] Ashburner | Nonlinear spatial normalization using basis functions[END_REF]. In most data analysis software packages, the reference image is the average T1 image provided by the Montreal Neurological Institute (MNI), which matches approximately the Talairach coordinate system [START_REF] Talairach | Co-Planar Stereotaxic Atlas of the Human Brain. 3-Dimensional Proportional System : An Approach to Cerebral Imaging[END_REF].

Voxel-based inference schemes are explicitly based on the assumptions that i) the functional images are properly co-registered, so that a location in the common space corresponds to the same region in the brain of each subject; ii) at a given spatial location in the reference space, the signal is normally distributed across subjects, so that the RFX and MFX statistics are Studentdistributed under the null hypothesis that no activation occurs. Both hypotheses might be wrong: the signal can be inhomogeneous across subjects [START_REF] Wei | Functional MRI of auditory verbal working memory: long-term reproducibility analysis[END_REF], so that normality assumptions are not met [START_REF] Thirion | Analysis of a large fMRI cohort: Statistical and methodological issues for group analyses[END_REF], and mis-registrations remain after spatial normalization of the datasets. The magnitude of such local shifts is probably 1cm in many brain regions (this can be observed for functional regions like the the motor cortex or the visual areas [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF], [START_REF] Stiers | Mapping multiple visual areas in the human brain with a short fMRI sequence[END_REF] or the position of anatomical landmarks [START_REF] Collins | Non-linear cerebral registration with sulcal constraints[END_REF]- [START_REF] Mangin | Object-based morphometry of the cerebral cortex[END_REF]). In addition, the number of subjects included in the analysis is generally small, so that RFX analyses are known to have a weak sensitivity.

In order to deal with the spatial mis-registration issues, most neuroscientists are thus accustomed to smoothing their datasets (8-12mm FWHM typically in group studies) to increase signal spatial overlap across subjects. This leads to biased and less precise localization of activated regions and may in some cases reduce sensitivity.

The interpretation of the boundaries of supra-threshold [START_REF] Mériaux | Combined permutation test and mixed-effect model for group average analysis in fMRI[END_REF] regions in group studies is not clear. Another approach consists in computing local or global anatomical warps that improve inter-subject co-registration [START_REF] Fischl | Cortical surface-based analysis II: inflation, flattening and a surface-based coordinate system[END_REF]- [START_REF] Wang | Brain surface parameterization using riemann surface structure[END_REF]. But such warps may require the additional use of anatomical landmarks, and it is not clear that different brains can be correctly warped onto each other. In particular, the variability in the large scale sulco-gyral anatomy [START_REF] Brett | The problem of functional localization in the human brain[END_REF], [START_REF] Rivière | Automatic recognition of cortical sulci of the human brain using a congregation of neural networks[END_REF] might imply that no such correspondences exist.

Note also that Talairach atlas was designed for subcortical structures.

In order to cope with non-normality of the signal across subjects, robust inference schemes, based e.g. on the sign test or Wilcoxon signed rank's statistic [START_REF] Mériaux | Robust statistics for nonparametric group analysis in fMRI[END_REF] have been designed. Moreover, permutation-based assessment of the group signal statistics [START_REF] Hayasaka | Validating Cluster Size Inference: Random Field and Permutation Methods[END_REF], [START_REF] Mériaux | Combined permutation test and mixed-effect model for group average analysis in fMRI[END_REF] yields an unbiased significance for the statistical maps across subjects, and thus bypasses some approximations implied by the use of random field theory [START_REF]Human Brain Function[END_REF].

However, performing a test on each and every voxel has a statistical cost (multiple comparison correction of the p-values), while many of these voxels are probably of little relevance to the cognitive function under study.

An interesting alternative is thus to perform inference at a higher level than the voxel level. In other words, one can consider functional regions, or structures, that are found active across subjects rather than active voxels. This point of view has been advocated by many groups that would use functional localizer paradigms to define brain regions before testing the activity of these regions in other conditions [START_REF] Saxe | Divide and conquer: A defense of functional localizers[END_REF]. In particular, regions of interest are frequently defined anatomically in order to ease functional studies [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] Eickhoff | Testing anatomically specified hypotheses in functional imaging using cytoarchitectonic maps[END_REF]. However, such regions are defined within a reference space (e.g. MNI space), which raises the aforementioned issue of misregistrations; moreover, such approaches define regions very coarsely [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] Eickhoff | Testing anatomically specified hypotheses in functional imaging using cytoarchitectonic maps[END_REF] (less than hundred regions for the entire brain). It is thus necessary to propose datadriven approaches.

In the literature, there is no generally accepted generative model of brain activity that could drive group inference procedures. Although few attempts have been

proposed recently [START_REF] Penny | Mixtures of general linear models for functional neuroimaging[END_REF]- [START_REF] Kim | Hierarchical dirichlet processes with random effects[END_REF], such approaches are likely to be confounded by the complexity of the data, the unknown extent and nature of the activations networks and the global cross-subjects variability. Therefore, a more pragmatic solution consists in modelling some structures of interest observed in the groups of subjects, and then to compare them in order to infer a grouplevel template of the observed data. Such approach are rule-based rather than based on a generative model of the data. Hereafter, such approaches will be called structural.

Structural approaches have to address several important questions:

• What are the structures of interest in each subject?

In the case of fMRI data, it is clear that the information of interest is coded in the maxima of activity maps, e.g. large supra-threshold clusters [START_REF] Hayasaka | Validating Cluster Size Inference: Random Field and Permutation Methods[END_REF], [START_REF] Simon | Automatized clustering and functional geometry of human parietofrontal networks for language, space, and number[END_REF],

scale-space blobs [START_REF] Coulon | Structural group analysis of functional activation maps[END_REF] or activity peaks [START_REF] Thirion | Finding landmarks in the functional brain: Detection and use for group characterization[END_REF]. Alternatively, some alternative approaches start with the prior definition of regions (parcels), based on clustering of anatomical and/or functional datasets [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF]. Some of these approaches might be somewhat coarse for a fine description of activated areas [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF],

[29], [START_REF] Thirion | Finding landmarks in the functional brain: Detection and use for group characterization[END_REF]. In this work, we rely on watersheds of supra-threshold areas, which is an intuitive and classical technique in pattern recognition [START_REF] Vincent | Exact euclidean distance function by chain propagation[END_REF].

• How to associate such regions across subjects ? This point may be more difficult, in particular because there exists clearly no isomorphism between individual active regions. While the position in a common space is an important information [START_REF] Simon | Automatized clustering and functional geometry of human parietofrontal networks for language, space, and number[END_REF], [START_REF] Thirion | Finding landmarks in the functional brain: Detection and use for group characterization[END_REF], [START_REF] Mériaux | Combined permutation test and mixed-effect model for group average analysis in fMRI[END_REF] there might be local variations that induce some ambiguities. In such cases, the relative position of neighboring regions might be of great importance [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF], [START_REF] Coulon | Structural group analysis of functional activation maps[END_REF]. In this work, we propose a relatively simple scheme to take this information into account.

• How to validate the sets of regions that have been associated across subjects ? In [START_REF] Coulon | Structural group analysis of functional activation maps[END_REF] a procedure that takes into account the individual feature quality, structural similarity between features and association strength, has been proposed. However, its complexity may be quite problematic for interpretation purposes. Here, we prefer to perform a spatial density test on the candidate regions, which allows a strict control on the specificity (type I error rate) of the method.

Finally, another important point is that structural methods involve many parameters in the modelling steps, and it is thus quite important to control the robustness of the results with respect to mild variations in the parameter setting.

In the present paper, we propose a framework that solves the aforementioned issues sequentially; in brief gives some means to account for and characterize intersubject differences, a key issue in group studies [START_REF] Thirion | Analysis of a large fMRI cohort: Statistical and methodological issues for group analyses[END_REF], [START_REF] Kherif | Group analysis in functional neuroimaging: selecting subjects using similarity measures[END_REF].

We describe the method in Section II, and some artificial and real benchmark datasets in Section III.

Importantly, our approach allows for an explicit control on specificity, which is shown in Section IV; in Section V we illustrate the improvement in terms of sensitivity and reliability of fMRI group analyses. Reliability is assessed by jackknife subsampling in a population of 102 subjects, and we show that the results of the proposed method are less dependent on the particular subgroup of subjects under study than standard voxel-based tests.

Finally, we describe the results of the method when applied to the whole group of 102 subjects. Technical issues and implications for neuroimaging studies are discussed in Section VI.

II. METHODS

A. Notations

Let us assume that a group of S subjects take part in an fMRI acquisition protocol while they undergo a certain cognitive experiment. After some standard preprocessing (distortion correction, correction of differences in slice timing, motion correction, normalization), the dataset of each subject is analysed in the General Linear Model (GLM) framework: for a given subject s ∈ {1, .., S}, let Y s be the dataset written as matrix (scans×voxel), and let X be the design matrix that describes effects of interest and confounds; the GLM proceeds by estimating the effect vectors β s such that

Y s = Xβ s + ǫ s , ∀s ∈ {1, .., S}, (1)
where ǫ s represents the residual matrix. The estimation is based on a maximum likelihood approach performed in each voxel, where the noise is assumed to be an AR [START_REF]Human Brain Function[END_REF] process [START_REF]Human Brain Function[END_REF], [START_REF] Worsley | A general statistical analysis for fMRI data[END_REF], [START_REF] Woolrich | Temporal autocorrelation in univariate linear modelling of fMRI data[END_REF]. Let c be the linear combination of [START_REF] Mériaux | Combined permutation test and mixed-effect model for group average analysis in fMRI[END_REF] the experimental conditions that is of particular interest;

c is also called a functional contrast. A certain statistic φ s can be computed in each subject s to assess the presence of a positive effect c T β s > 0 in each voxel of the dataset, e.g.

φ s (v) = E(c T β s (v)|Y s ) var(c T β s (v)|Y s ) (2)
at each voxel v.

Our method takes as input the activations maps φ s of each subject s ∈ {1, .., S}, which can be thresholded at a certain significance level, using either voxel-level or cluster-level assessment. In what follows, we assume that, given the significance level P, there exists a known threshold θ 0 such that P (φ s (v) > θ 0 |H 0 ) < P at any voxel v, where H 0 represents the null hypothesis that no activation is present. Our analysis procedure is illustrated in Fig. 1 and consists of four steps, which are detailed in the next parts.

[Figure 1 about here.]

B. Intra-subject Structural analysis

In the absence of a sound prior on the nature or the position of activated foci, our approach first extracts regions of interest in each dataset. It seems particularly meaningful to segment the main peaks of activity within the supra-threshold components of the statistical maps (φ s ) s∈{1,..,S} : the connected supra-threshold components in each subject s are thus segmented into I(s) regions using a watershed method, so that each segmented region is associated with a local maximum of the map φ s . Let (a s i ) i=1..I(s) be the corresponding maxima for subject s, and (t s i ) their MNI coordinates (which approximate Talairach coordinates). It is useful to have a graphical representation of the spatial relationships between the segmented regions in each subject (see Sec. II-D). Several models may be used to build such a graph, for instance the neighboring relationships between adjacent regions. The resulting graphs (G s ), s ∈ {1, .., S} are undirected, and they may contain cycles.

We propose an alternative, that produces acyclic graphs:

The list of maxima of any connected regions can be organized according to the order relation O: a s i Oa s j if and only if the corresponding regions are neighboring and if a s j is the highest maximum in the vicinity of a s i . The set of these structures across regions defines a directed acyclic, possibly disconnected graph in each subject. We consider the undirected graphs with the same edges; these will be denoted G s , s ∈ {1, .., S}.

Assuming that the activated regions are aligned along some sulci, hence in one-dimensional structures, the treelike representation given by G s may code quite well their spatial organization.

Note that the normalization procedure can, and in principle should, take place after the intra-subject analysis part. We have implemented both solutions on a real dataset, and have not noticed any significant difference in the global outcome of the method.

C. Spatial statistics

Given that the initial threshold P should preferably be kept low to avoid false negatives, the first step necessarily results in several false positives. A statistical test on the spatial distribution of the maxima is thus performed to control the false positive rate. Only regions with across-subject reproducible activity are of interest.

Thus we build a spatial statistic to remove the local maxima of each subject that are far from local maxima of other subjects. This spatial statistic is the density of presence of supra-threshold local maxima in the other subjects. Let τ = (x, y, z) be a position in the common [START_REF] Mériaux | Combined permutation test and mixed-effect model for group average analysis in fMRI[END_REF] space. For each subject s ∈ {1, .., S}, we define

D s (τ ) = σ∈{1,..,S}-{s} I(σ) i=1 exp - τ -t σ i 2 2δ 2 τ (3)
The parameter δ τ represents an inter-subject spatial variability and is set to 10mm.

The distribution of the quantity D s (τ ) can then be compared, in every location, to its distribution under the null hypothesis H 0 . The null hypothesis that we consider here is global, i.e. it means that there is no specifically task-related region in the brain. Under this assumption, the spatial density of local maxima is uniform in the brain volume. We estimate the distribution of D s (τ )

under H 0 by random resampling of the position of the activation maxima (a σ i ) i=1..I(σ),σ =s within the brain volume. Let Ds be the surrogate distribution obtained after k resamplings (k = 10 typically).

Then, let α be a significance level, and let u α be a threshold on the values of D s such that P (D s (τ ) > u α |H 0 ) < α; an estimator of u α is given by the α-quantile of the density Ds

u α = arginf u 1 Ω I Ds>u (τ )dτ < α (4)
where Ω = dτ is the brain volume.

Importantly, the test is performed for a small number of spatial locations (t s i ) i=1..I(s) . Hence its significance can be corrected using a Bonferroni procedure, i.e., by replacing α by α I(s) in Eq. ( 4). An example is provided in Fig. 2.

One might be concerned with the behaviour of the method, assuming that the null hypothesis has been rejected in some regions of the brain: Does the test remain valid in the other regions, given that the global null hypothesis of a uniform density of maxima has been rejected ? In fact, in such case, the resampled distribution Ds is an overestimation of the true null distribution under the null hypothesis P (D s |H 0 ), which means that the ensuing test is conservative, hence valid. This fact is evident in Fig. 2, where the (null) mode of the resampled distribution is shifted to the right, with respect to the mode of the non-resampled distribution.

[Figure 2 about here.]

The test is iterated in all the subjects, then nonsignificant maxima at the desired significance level are rejected. The process can be iterated in order to refine the spatial model. Let I(s) ≤ I(s), s ∈ {1, .., S} be the number of remaining regions in each subject. Since we control the probability that one false positive region might show up in a given subject at level α, given ν ∈ {1..S}, the probability of one false positive region in ν subjects over S is controlled by the binomial law B(ν, S, α).

D. Finding correspondences using Belief Propagation networks

The statistical procedure leaves us with a set of candidate regions which are spatially clustered across subjects. Then the core part of the procedure consists in finding which regions correspond across subjects. This problem is probably the most difficult one, since one would like to obtain explicit correspondences between individually segmented regions, although there cannot be a one-to-one correspondence across subjects. For instance, a couple of neighboring active foci in one subject might not be distinguishable in some other subject. Our solution consists in estimating, for each pair (s 1 , s 2 ) of subjects, the probability that the region a s1 j in subject s 1 is the analogue of region a s2 i in subject s 2 . Then, given these probabilities, cliques of cross-subjects regions will be constructed; this will be detailed in section II-E.

We search for inter-subject correspondences in the relative positions of activated areas. While this does May 29, 2007 DRAFT not rely on a physical model, it builds on the heuristic that across subjects, positions of activated regions should be locally similar, though not identical in the common space. For instance, it is logical to favor configurations in which couples of neighboring regions with similar relative positions in two subjects will be grouped in two different cliques according to their relative position.

We base our search on a graphical model of the position of the maxima in each subject. This model is provided either by the undirected acyclic graph G s or the undirected and possibly cyclic graph G s defined in section II-B, from which non-significant nodes have been removed. Probabilistic associations are then searched for each pair of subjects, using a belief propagation (BP) algorithm [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]. Given a reference subject s 1 and a target subject s 2 , the associations are initialized as

P (a s2 i ← a s1 j ) ∝ exp - t s2 i -t s1 j 2 2δ 2 τ (5)
with appropriate normalization, where P (a s2 i ← a s1 j ) stands for the probability that the maximum a s2 i in subject s 2 is the homologue of maximum a s1 j in subject s 1 . These probabilities are refined by belief propagation; for each edge (jk) of the graph G s1 or G s1 , messages are sent from a s1 j to a s1 k to quantify the probability of association between a s1 k and (a s2 i ) i=1..I(s2) :

m jk (i) ∝ I(s2) l=1 P (a s2 l ← a s1 j ) exp - (t s2 i -t s2 l ) -(t s1 k -t s1 j ) 2 2δ 2 τ (6)
with appropriate normalization (

I(s2) i=1 m jk (i) = 1).
Eq. ( 6) simply means that whenever the positions of local maxima t s2 l and t s2 i in subject s 2 and t s1 j and t s1 k in subject s 1 form a parallelogram, the configurations are favored in which a s2 l and a s1 j on the one hand, a s2 i and a s1 k on the other hand, are associated. As shown in Fig. 3, taking into account the relative positions of the maxima improves cross-subjects correspondences by compensating global translation effects.

[Figure 3 about here.]

The beliefs P (a s2 i ← a s1 j ) and messages are then updated and normalized according to the formal laws of BP [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF], [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF]. Note that the graphs G s have no loops, so that convergence is straightforward. We have also used loopy belief propagation, based on the graphs G s , which did not raise any issue concerning the convergence of the correspondence probabilities. Furthermore, the choice of G s or G s was not found to be crucial in the method. In our experiments, we use G s by default.

The estimation of the probabilities is performed on each pair of subjects in the group. As an important note, all the quantities used here are asymmetric. In particular, the graphs G s1 and G s2 have a priori different structures, so that the probabilities P (a s2 i ← a s1 j ) and P (a s1 j ← a s2 i ) might be quite different after convergence. This is particularly obvious in the case of many-to-one correspondences, given that the number of maxima in a given region may vary a lot across subjects. The next step essentially chooses which of these correspondence are meaningful at the population level.

E. Extracting homologous regions

All the probabilities of all pairwise associations between subjects are then arranged in a common belief matrix B. A row of B contains all the probabilities that a maximum a s i of a subject is associated with all maxima (a σ j ) j=1..I(σ),σ∈{1,..,S} of the other subjects; note that the associations with the other maxima (a s j )j = 1..I(s) within the subject s itself are null. We can also interpret B as the adjacency matrix of a probabilistic association graph between regions that are pooled across subjects. May 29, 2007 DRAFT Some of these associations may not be realistic, but in that case the two maxima should not have strong common associations with other maxima. To deal with such cases, we proceed with the extraction of the maximal cliques of the belief matrix, i.e. groups of maxima that have mutually strong associations. In our case, the association probabilities are asymmetric, thus the maximal clique approach requires that association probabilities are relatively high bidirectionally, in a sense detailed thereafter.

Many clustering procedures are possible, e.g. hierarchical clustering techniques, using average or maximal linkage heuristics or replicator dynamics (RD). Since the latter procedure is more data-driven (it does not require a prior definition of the number q of clusters to be found), we describe it in more details, but we also suggest to use average-link agglomerative clustering, where the number of desired clusters is q = mean s I(s) the average number of regions per subject. Note that in that case, the association probabilities are symmetrized.

A formal definition of graph-theoretical cliques is

given in [START_REF] Pavan | Dominant sets and pairwise clustering[END_REF], in the case where the matrix B is sym- neighbors. Finally, it is shown in [START_REF] Pavan | Dominant sets and pairwise clustering[END_REF] that it amounts to define a membership vector x on the graph vertices, and then to solve the program maximize x ′ Bx subject to x ≥ 0 and x ′ u = 1 [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF] where u is the vector of ones with the same size of

x. Finally, still in the case where B is symmetric, this problem can be solved using replicator dynamics equations to B (see e.g. [START_REF] Lohmann | Using replicator dynamics for analyzing fMRI data of the human brain[END_REF], [START_REF] Neumann | Meta-analysis of functional imaging data using replicator dynamics[END_REF]). Replicator dynamics consist in initializing, randomly or not, a positive vector x (0) whose length is equal to the total number of vertices, and then in iterating the update rule

x (i+1) = (Bx (i) ).x (i) x (i) ′ Bx (i) (8)
where * . * stands for the element-wise product. After a few iterations, almost all the components of x vanish, and the other ones correspond to a maximal clique of the belief graph. The clique is removed, and the process is repeated until no non-trivial clique is found. Other rules than the replicator dynamics can be used instead [START_REF] Pelillo | Payoff-monotonic game dynamics and the maximum clique problem[END_REF], but we experienced that Eq. ( 8) works efficiently.

We noticed that this procedure tends to over-segment the graph B, which is natural due to the restrictive definition of the maximal cliques (or dominant sets, see above), but this is not a problematic issue, since further merging of cliques remains possible. For instance, the procedure can be iterated based on cliques instead of regions, thus yielding larger cliques.

Finally, all the cliques that contain maxima from at least ν (e.g. S 2 ) subjects over S are retained.

F. Derivation of a group template

This procedure provides us with clusters of activated regions defined across subjects. It does not require that all subjects are represented for a given activated region, and therefore is able to account for some inter-individual differences. A large cluster means that subjects typically have an activated region that corresponds to this cluster.

In order to make group maps, we assume that the positions t i s of the maxima within each clique are normally distributed, and thus represent the cliques through their 95% confidence regions (CR) in the common (MNI) and random or mixed-effects analyses, which yield the probability that the effect observed in any subject of the population will be positive.

G. Parameters and implementation issues

The method requires few prior parameters: the initial threshold of activity maps P and especially the spatial pvalue α are chosen in order to control the number of false positives. We take typically P = 0.001, uncorrected for multiple comparisons, and α = 0.2, but α can be tuned to obtain explicit confidence levels depending on ν (see below). The other parameters are the spatial relaxation distance δ τ , which we set as a typical inter-subject variability magnitude δ τ = 10 mm, and the number ν of subjects required for the final selection of cliques. ν is important since it explicitly controls the reproducibility of a region across subjects. For instance ν = S 2 yields regions that can be expected to be found in half of the subjects. Note that all these parameters (P, α, δ τ , ν) might be changed reasonably without creating inconsis-tencies.

The control of false positive regions is based on the control performed in each subject in section II-C.

For instance, if one were controlling the rate of false detections, Eq. ( 4) provides us with the probability of one false alarm on any subject, so that an upper bound of the probability of forming a clique with regions from ν subjects over S under the null hypothesis is given by the binomial law B(ν, S, α): hence, the probability of getting at least one clique of ν subjects, among S is given by

p < n≥ν B(n, S, α) (9)
The computation time of the method, implemented with a C/Python code based on numpy and nipy environments (http://projects.scipy.org/neuroimaging/ni/), is about one minute for a dataset of 10 subjects. The complexity of the method is roughly quadratic, given that all pairs of subjects are submitted to the algorithm described in Section II-D. However, the main bottleneck with large datasets is the clique extraction procedure (section II-E), and the computation time is largely datadependent. For a dataset of 100 subjects, the proposed method take approximately 1h on 3GHz Pentium IV PC running Linux.

III. EVALUATION DATA

A. Synthetic Noise

The algorithm was tested on synthetic datasets with no simulated activations, in order to ensure that the false alarm rate was controlled using equation [START_REF] Collins | Non-linear cerebral registration with sulcal constraints[END_REF]. The simulations is reported, and compared to the theoretical bound provided by Eq. ( 9).

B. Artificial Activations in Synthetic Noise

We also applied the procedure on synthetic data with activation added to the correlated noise.

C. Real fMRI data

We used an event-related fMRI paradigm that comprised ten experimental conditions. Subjects were presented with a series of stimuli or were engaged in This pre-processing was performed using the SPM2 software (www.fil.ucl.ac.uk, [START_REF]Human Brain Function[END_REF]). In particular, spatial normalization was performed using default parameters (non-rigid, low frequency deformation with 8*8*7 basis functions [START_REF] Ashburner | Nonlinear spatial normalization using basis functions[END_REF]); the normalized images were checked in all the subjects to prevent any gross mistake in the image May 29, 2007 DRAFT co-registration. A slight smoothing was performed (5mm FWHM). Standard statistical analysis were also carried out with SPM2, using the usual high-pass filtering and AR(1) whitening.

In the study of the group data, we concentrated on a specific cognitive contrast that shows the activation elicited by the computation task, after video or audio instruction, from which the mere sentence reading/listening effect is subtracted. Such a cognitive contrast is assumed to yield areas specifically activated in the computation task.

In parallel, the grey/white matter interface was segmented in each subject using the Brainvisa software (http://brainvisa.info/), and is used for rendering.

D. Assessment of the reproducibility

Controlling the specificity of the analysis is not sufficient to have reliable brain maps; another concern is to control the risk of overfit in small populations that could result in poor generalization of the regions found to other groups of subjects.

We dealt with this concern by performing the above analysis in ten disjoint groups of 10 subjects sampled from a set of 102 subjects present in the database. We computed an inter-group reliability index by analyzing how often a voxel is declared jointly active across groups, using the procedures described in [START_REF] Thirion | Analysis of a large fMRI cohort: Statistical and methodological issues for group analyses[END_REF], [START_REF] Liou | Bridging functional MR images and scientific inference: Reproducibility maps[END_REF],

[42]: The reproducibility map that gives the number of times each voxel is declared active is computed, and its histogram is derived; this histogram is modelled by a mixture of two binomial distributions, and the index κ ∈ [0, 1] measures the accordance of the bimodal model with the data, which in turns reflects the coherence of the binary maps given as input to the model. If κ is close to 0, there is a very little agreement on which voxels are active, while there is a very good agreement if κ is close to 1.

The reliability was estimated from 100 different random splits of the group. This was computed for Random effects analysis (RFX), the same RFX analysis after 12mm FWHM smoothing of the data (SRFX), a Mixed Effects (MFX) analysis [START_REF] Worsley | A general statistical analysis for fMRI data[END_REF], an RFX analysis thresholded at the cluster level (CRFX), a Parcel-based RFX (PRFX), and our estimate of confidence regions (CR). More precisely, RFX, SRFX, PRFX and MFX maps were thresholded at the uncorrected p < 0.001 level; the CRFX map was built by taking the voxels with a signal significant at p < 0.01 uncorrected level, then clusters of connected supra-threshold voxels were formed, and further selected if their size was significant, at p < 0.05 corrected level [START_REF] Hayasaka | Validating Cluster Size Inference: Random Field and Permutation Methods[END_REF]. Lastly, the CR maps contain the 95% confidence regions for the presence of maxima in ν = 4 over 10 subjects, with α = 0.2. Note that the PRFX procedure is performed as in [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF], and parcels are recomputed for each randomization and sample. The parameters were chosen in order to guarantee that the specificity of the different methods is roughly equivalent, and that the parameters correspond to standard choices.

Additionally, we compared the results of the CR extraction procedure with different graphical models, G s , G s , or the trivial graph with no link (i.e. without the belief propagation algorithm). For this purpose, we used a functional contrast that showed region involved in the processing of auditory instructions, because this functional contrast elicits many neighboring activation foci (a s i ) (about ten in average) in each temporal lobe and each subject.

E. Analysis of the group data

We computed the RFX map and extracted Confidence Regions of the areas activated by the computation task across subjects. Then, we computed the average signal May 29, 2007 DRAFT and position of the regions in each subject, whenever they are defined. We performed some data-driven clustering of these profiles, which yields an assessment of the population homogeneity. Then we regressed this data against side information that was obtained form the subjects: in this case, we used the age of the subjects, their sex, and their ability to perform the mental rotation in 3D of an object, measured by the rate of correct response in a psychological test. The regression procedure reads simply:

Y g = X g β g + ε (10)
where Y g is a (subjects, voxel/ROIs) data matrix that represents the average ROI-based or the voxel-based activation signal, X g a group-level design matrix of size (subjects, regressors) that represents the covariates of interest across subjects, β g the second-level regression parameters of size (regressors, voxels/ROIs) and ε the residual. Then the voxel-or ROI-based significance of β can be assessed using a standard t-test.

Finally, in the case of the ROI analysis, not only the voxel-based average signal, but also the cross-subject ROI position in the common space can be given as input to model [START_REF] Hellier | Retrospective evaluation of intersubject brain registration[END_REF]. In this case, a chi square test can be used to assess the correlation of the regressors in X g with the ROI positions.

IV. SIMULATION RESULTS

A. Controlling the False Positive Rate in Synthetic Noise

We have applied the procedure to n=100 synthetic datasets generated as detailed in section III-A, with different values of the parameter α. The detection rate is compared with the theoretical bound given in Eq. ( 9).

The results are reported in Fig. 4, which shows that the control of false positives is conservative, since the rate of detected regions obtained is below the predicted value.

In fact, the control is probably too conservative. The only case where the control may be problematic is that for very low values of α, the correct definition of u α (see Eq. ( 4)) requires a very accurate estimate of the right tail of D s under the null hypothesis, hence many resamplings. In practical cases, we found 10 resamplings to be sufficient.

We have repeated the procedure with different values of the first-level threshold P, from 0.05 to 10 -4 .

Although the number of false positive depended on the particular value of P that was chosen, it always remained under the theoretical bound (which does not depend on P).

[Figure 4 about here.]

B. ROC Analysis in Synthetic Activation and Noise

We computed the ROC curve on synthetic dataset with embedded activation in four regions. The exact position of each activated region may vary from 0 to 3.5 voxels, while the amplitude of the response is allowed to vary around a mean of -6dB or -10dB. ROC curves, that represent specificity/sensitivity compromise when detection parameters vary, are presented in Figure 5 for RFX analysis in the proposed CR method, in the different situations. Note that we consider only the part of the curve with a specificity control below 0.01, since weaker controls are of no practical interest. The simulations are repeated 100 times.

[Figure 5 about here.]

Except for the situation with no jitter, the CR method clearly outperforms the RFX method. In addition, the CR method is clearly less sensitive to the SNR level than the RFX statistic. In the no jitter case, the CR method performs better at high specificity levels, but reaches a [START_REF] Mériaux | Combined permutation test and mixed-effect model for group average analysis in fMRI[END_REF] plateau at a lower level when the specificity control is weaker. The reason is that the CR is not meant to achieve 100% sensitivity at the voxel level, since the extent of the CR corresponds to the inter-subject variability of the areas. By contrast, the RFX statistic actually tests the presence of activation in each voxel, and thus can detect 100% of the activated voxels in the absence of jitter.

V. RESULTS ON REAL DATA

A. Gain of sensitivity in a small group of subjects

First, we performed a group analysis in a small group of 10 randomly chosen subjects within the whole group.

A CR map, obtained at a p < 0.05 significance level for the computation-specific contrast is presented in Fig. 6, together with an RFX map, thresholded at a p < 10 -3 , uncorrected, or at the cluster level at p < 0.05, corrected. When the RFX map is thresholded at the same significance level (p < 0.05, corrected for multiple comparisons) at the voxel level, no voxel survives the thresholding procedure.

[Figure 6 about here.]

In spite of the strong type 1 error control, the CR map contains 19 significantly active regions. In particular, it clearly shows symmetric parietal regions involved in the computation task, while these regions are not detected with the RFX procedure (voxel-or cluster-level statistics). This is in agreement with the literature [START_REF] Lemer | Approximate quantities and exact number words: dissociable systems[END_REF] (see also Fig. 9).

In order to test the robustness of the method, we have repeated the experiment with unsmoothed/smoothed data (FWHM=10mm), with higher or lower first-level thresholds P = 10 -2 , 10 -3 , 10 -4 , with the spatial normalization before or after the watershed, or with a small spatial jitter of activation images across subjects (one voxel). All these changes had a very weak impact on the resulting CR maps. For instance, to test the robustness to small spatial shifts, we randomly shifted the datasets from zero or one voxel in one direction (x,y or z) with equal probability (1/7). We performed the RFX and CR analysis of the resulting group data and derived the reproducibility index κ (see Sec. III-D) from ten such group maps. Over 100 repetitions, we obtained a mean value of κ = 0.61 (range 0.58 -0.64), which has to be compared with a mean value κ = 0.42 (range 0.37 -0.47) in the case of the RFX map.

Importantly, the CR procedure not only provides a group-level activity map, but-also explicit correspondences between active regions at the subject-level and the group data. This is illustrated in Fig 7.

[Figure 7 about here.]

B. Between-Group Reproducibility

The reliability of the group analysis method, assessed by reproducibility index estimated across voxels for different splits of the populations into groups of 10 subjects is shown in Fig. 8. This shows that the proposed CR method outperforms RFX, MFX and SRFX, and to a lesser extent, PRFX and CRFX. It is important to note that this procedure is based on the active or inactive status of each voxel, and should not be favorable to non voxel-based analyses a priori. Moreover, the present results are not related to the particular choice of thresholds or p-values, and other reproducibility indexes also yield similar effects (not shown).

[Figure 8 about here.]

Moreover, for a contrast that shows regions involved in the processing of auditory instructions, we found that the CR regions were more reproducible when obtained using either the acyclic graph G s (κ = 0.589 ± 0.01) [START_REF] Mériaux | Combined permutation test and mixed-effect model for group average analysis in fMRI[END_REF] or the spatial connectivity graph G s (κ = 0.597 ± 0.01) than a trivial graph without edges (i.e. without the belief propagation algorithm, see Eq. ( 6), κ = 0.539 ± 0.01).

C. Analysis of a large population

The application of the CR methods to the whole group of 102 subjects yielded q = 45 regions, with a corrected p-value of 0.05 of making one false detection. The CR map is presented with the RFX map on these same subjects in Fig. 9. As is usual with large sample sizes, the RFX map shows very wide activated areas. This is is the result of blurring process inherent to the inter-subject variability.

[Figure 9 about here.]

We describe the properties of the 45 resulting clusters in Table I. In particular the anatomical labels of the regions are found in [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].

[Table 1 about Next, we computed the average signal per region per subject, and tried to characterize the population by unsupervised classification techniques. We call the average signal per region for each subject the profile.

Based on simple Euclidean distance between profiles, we have performed some agglomerative clustering of the population, using an average linkage procedure. The results are shown as a dendrogram of the subject's profiles in Fig. 10. In this case, it clearly shows that the population is mainly divided into one group of 97 subjects, and 5 isolated subjects. One can conclude that the population is rather unimodal, with a few outliers. A closer inspection of the outlier datasets reveals that four of them had no significant activations, and the last one had an odd pattern of activity, probably confounded by motion or another low-level artifact.

[Figure 10 about here.]

Finally, we regressed the voxel-based activity maps as well as the profiles against three regressors of interest defined in each subject: age, sex, and ability to perform a 3D mental rotation (see Eq. 10). We found no effect of age in either case 1 .

Concerning sex, we found in the voxel-based analysis a region where the magnitude of the activity is larger for males than for females [z = 5.35, p < 0.05, corrected at (18, -68, 60)mm]. As this place is on the posterior edge of the parietal lobe, and not in a significantly activated region, the interpretation of this result is quite unclear. The ROI-based analysis revealed that there was indeed a positive effect on one ROI [z = 2.73, p < 0.01, uncorrected at (15, -68, 52)mm], but moreover that there was a significant effect of the sex on the ROI coordinates in MNI space across subjects [χ 2 3 = 13.95, p < 0.01, uncorrected at (15, -68, 52)mm], indicating that there might be some systematic shift effect between males and females. Importantly, no such effect can be observed using the voxel-based analysis. [Figure 11 about here.] In summary, while voxel-based and ROI-based analyses show similar effect (in the statistical sense) of the age, sex or 3D task performance on the fMRI signal, the ROI-based analysis clearly indicates that these differences could be spatial, thus possibly anatomical, rather than merely quantitative.

VI. DISCUSSION

Incrementing our knowledge on human brain function requires the common analysis of datasets from different subjects. While standard analyses assess the significance of effects at the voxel level, we show here that this procedure is not optimal since it suffers from the heterogeneity of the signal measured in different subjects, and from mis-registrations. To deal with these issues, taking into account the absence of a satisfactory generative model of brain activity in groups of subjects, we presented a rule-based, structural approach that extracts structures of interest in each subject's dataset and builds a group model from the structures of each subject.

Our solution is in the same spirit as a previous structural approach [START_REF] Coulon | Structural group analysis of functional activation maps[END_REF] based on the detection of scale-space blobs and the discovery of correspondences with a Markov Random Field (MRF). Our definition of activated regions by watershed analysis of supra-threshold regions is simpler and spatially better defined than the scale-space blobs. Our Belief Propagation scheme is quite comparable while simpler (see Fig. 3) than the MRF model [START_REF] Coulon | Structural group analysis of functional activation maps[END_REF], since the latter had to take into account some idiosyncrasies of the scale-space blob model. Furthermore our procedure inherits the good convergence properties of BP algorithms [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF]. For the representation of the spatial structure of activated regions, we noticed that a loopy BP algorithm performed as well, or even slightly better than a tree-based BP algorithm; in any case, introducing the BP scheme markedly improved the reliability of the correspondences across groups with respect to a standard approach based only on the position in the common space (i.e. on Eq. ( 5)).

Other alternatives to our procedure are the clusterbased inference (CRFX, [START_REF] Hayasaka | Validating Cluster Size Inference: Random Field and Permutation Methods[END_REF]) and parcel-based inference (PRFX, [START_REF] Thirion | Dealing with the shortcomings of spatial normalization: Multi-subject parcellation of fMRI datasets[END_REF]). In particular, it is shown in Fig. 8 that these two are almost as reliable as the described approach. They suffer, however, from important drawbacks: cluster-based inference assumes that only wide suprathreshold clusters are worth reporting, which is not always true (see Fig. 6). Parcel-based inference that consists in making R/MFX tests on parcels instead of voxels, builds parcels of arbitrary size, and thus does not always correctly model the fine-scale activation pattern in each subject. The present approach based on watershed analysis of activity maps might reveal finer scale details (see Fig. 7). This might become especially important with the advent of high-resolution fMRI acquisition techniques [START_REF] Grill-Spector | High-resolution imaging reveals highly selective nonface clusters in the fusiform face area[END_REF]. Compared with the parcel-based approach for which the inter-subject correspondence is assumed a priori, the described solution relies on an a posteriori scheme that yields more information on the May 29, 2007 DRAFT reproducibility of a spatial pattern.

The specificity control of our procedure is correct on surrogate data (see Fig. 4), though much too conservative in many instances. This is due to simplifying (and conservative) assumptions used for the derivation of the test (see Sec. II-G). Finding a tighter upper bound of the error rate might be an important topic for future studies. Although first level statistics are used to define the regions of interest across subjects, an important point is that the control of false positive regions at the group level (see Eq. ( 9)) does not explicitly take the first-level statistics into account. In fact, the whole procedure, and in particular the statistical test described in section II-C is a region selection and association procedure that is blind to the first-level definition procedure (here, a p-value thresholding and watershed separation), and remains valid as along the first-level procedure is performed independently in each subject. In particular, it adapts quite automatically to the variability and the noise level in the dataset. Moreover, it should be stressed that first level statistics are usually not very reliable, because unmodeled effects (physiology/motion) can have a great impact on the effect significance estimation.

Our solution is well adapted to the characteristics of fMRI data, in the sense that it optimizes the compromise between sensitivity and specificity, as shown by the ROC curves in Fig. 5. The CR method outperforms the RFX thresholding procedure in most instances, especially when the homogeneity between subjects is low in terms of spatial or quantitative functional information. The only case where RFX outperforms the proposed method is in the absence of jitter, and for a weak control of specificity: in practice, both assumptions are unrealistic.

For the analysis of a standard group (10-15 subjects), see Fig. E) tends to oversegment the active regions, but this effect can be solved using hierarchical clustering instead. This problem might also be by-passed in the future e.g. by using multi-scale methods [START_REF] Coulon | Structural group analysis of functional activation maps[END_REF], [START_REF] Poline | Analysis of individual brain activation maps using hierarchical description and multiscale detection[END_REF]. It is important to note however that this effect is not significant for small sample sizes.

Our aim was also to refine the conventional point of view on the localization problem in fMRI data analysis [START_REF] Brett | The problem of functional localization in the human brain[END_REF]. In particular, an activity map of one subject in its native anatomical space is not comparable to a group map presented on an average anatomy: in the latter case, we present the locations where individual subjects drawn from the group typically activate, while in the former we present for a specific subject the regions with significant activity (see Fig. 7). This should help to better interpret fMRI group studies results. In standard analyses, it is impossible to distinguish between regions for which all subjects show a small increase of activity from regions for which only some subjects demonstrate increased activity.

Moreover, finding correspondences across subjects allows us to make statements on the dissimilarity between subjects [START_REF] Kherif | Group analysis in functional neuroimaging: selecting subjects using similarity measures[END_REF], which is another blind spot of traditional M/RFX studies. Characterizing inter-subject differences in an interpretable way is essential if neuroimaging data is to be compared with genetic or behavioral information.

An example is given in Sec. V-C, where one can make some inference on the between-subject variability by trying to explain differences in size/and or position of regions across subjects by some information that is available on these subjects. This kind of inference is possible, but quite cumbersome in the traditional voxel-based domain, due to a curse of dimensionality (the number of voxels is too high), and because it is not clear whether the variability can be attributed to differences in the signal level across subjects or to the position of the regions. Our procedure, for instance, indicated that the positions of the active regions, thus the functional anatomy, plays a non-negligible role in group discrimination. This point will be further studied in the future using e.g. multivariate classification/regression analysis techniques.

Another important question is whether this kind of analysis can be generalized to group comparison, which is important e.g. for the characterization of brain diseases. Although the answer is probably case-dependent, one possibility consists in pooling the subjects to define the ROIs, then derive subject or group profiles as in Sec. V-C, and to study possible group differences at the ROI level. However, we acknowledge that some cases may be problematic, e.g. if the inter-group differences reduce the sensitivity of the spatial test to discover some of the ROIs. On the other hand, the proposed method provides the opportunity to compare the regions position or shape across groups, which is not afforded by voxelbased models.

The present work is also an attempt to automatically find correspondences across subjects by associating activated areas with close relative positions across subjects.

In particular, watershed analysis of the individual maps is used to define target regions. In order to enhance the understanding and interpretation of inter-subject variability, future developments might consider the use of generative models in the spatial and/or signal domain, based e.g.

on Dirichlet Process Mixture Models of the fMRI data [START_REF] Kim | A nonparametric bayesian approach to detecting spatial activation patterns in fmri data[END_REF], [START_REF] Kim | Hierarchical dirichlet processes with random effects[END_REF]. Alternatively, a more anatomical point of view may be introduced, e.g. by defining the position of ROIs with respect to macro-anatomical features (sulco-gyral May 29, 2007 DRAFT anatomy [START_REF] Rivière | Automatic recognition of cortical sulci of the human brain using a congregation of neural networks[END_REF]). This opens the way to an object-oriented representation of the functional anatomy [START_REF] Mangin | Object-based morphometry of the cerebral cortex[END_REF]. It is also worthwhile to note that MNI/Talairach space does not play any particular role in the present method, so that any valid -and non-Euclidean -normalized space such as the one on the cortical surface [START_REF] Fischl | Cortical surface-based analysis II: inflation, flattening and a surface-based coordinate system[END_REF], [START_REF] Clouchoux | Anatomically constrained surface parameterization for cortical localization[END_REF] can play the same role.

VII. CONCLUSION

In this work, we have shown that describing and comparing datasets with high level information instead of the usual voxel-based activity may benefit both the sensitivity and the reliability of fMRI group analyses.

Moreover, this approach does not lose the spatial information through an averaging process, but enables neuroscientists to make explicit comparisons between subjects or groups of subjects in a more rigorous conceptual setting.

Fig. 1. Flowchart of our method for structural analysis of group data. This is a pipeline, or set of procedures, that produces a group-level representation of individual activation patterns based on reproducibility analysis; we illustrate it for a motor activation study. The input to the method consists of activation images, one for each subject. An intra-subject structural analysis is first performed, resulting in a set of activated regions. The cross-subject spatial density of activated regions is derived, and only the maxima that fall in the highest density regions are further considered. Probabilistic correspondences are then found between the regions of each pair of subjects, using a belief propagation algorithm. Finally, associated regions are segregated into inter-subjects cliques, so that each region of reproducible activity is labeled consistently across subjects, and confidence regions for the position of these ROIs are derived at the group level. In that case, t s 1 i = 0, t s 1 j = 1 and t s 1 k = 2, while t s 2 l = 0.7, t s 2 m = 1.7 and t s 2 n = 2.7; δτ = 1.4. The related graphs G s 1 and G s 2 are shown below; in this case they are isomorphic. The associations are initialized using Eq. ( 5), and then refined using Eq. ( 6): Clearly the message passing algorithm enhances the probabilities P (a s 2 l ← a s 1 i ), P (a s 2 m ← a s 1 j ), and P (a s 2 n ← a s 1 l ), thus compensates the effect of the global translation between the two datasets. This effect is also present with the converse probabilities P (a s 1 . ← a s 2 . ). Fig. 5. ROC curves for the RFX and CR maps for different SNR levels and across-subject jitter magnitudes. (a) with no jitter of the activation position, ROC curves are presented for the RFX (blue) and CR method (red), for a mean SNR of -6dB (continuous line) or -10dB (dashed line). (b) and (c): the same curves, with a jitter of magnitude of 1 voxel in each direction (b), or 2 voxels in each direction (c). Except for the situation with no jitter, the CR method outperforms the RFX method. In the case of no jitter, the CR method performs better at high specificity levels, but plateaus at a lower level when the specificity control is weaker.

  metric. They are termed the dominant sets of the graph, and their definition relies on two conditions: i) that the similarity value of each element of the clique should be high enough with respect to the average similarity of the other elements ii) that any element outside the clique should have a weaker similarity with the clique than the elements of the clique, where the similarity values are computed from the affinity values between graph

  noise only datasets are meant to simulate multi-subject activation maps in which no specific region is activated. Masks of the brain volume were extracted from S = 10 subjects in the true dataset; the maps were filled with random normally distributed values and slightly smoothed May 29, 2007 DRAFT (FWHM = 1.17 voxels, i.e. 3.5mm at 3 mm resolution) in order to mimic the intrinsic spatial correlation of fMRI data. The signal magnitude was then corrected to have variance 1 in each dataset. The whole procedure was applied to these maps for values of the parameter α ranging from 0.1 to 0.8 in 0.05 steps, with δ τ = 10mm, ν = 5 and P = 10 -3 . We also tried with different values of P, ranging from 0.05 to 10 -4 . Finally, the number of detected regions in the group of subjects, over n = 100

  Four distant regions are added with some signal in order to model spatially coherent activity in the group of subjects. The size of the activated regions varied from 20 to 50 voxels. According to the simulation, their mean position (center of mass) was jittered with a magnitude of 0, 1.7 or 3.4 voxels standard deviation, which represents 0 to 10.4mm at 3mm resolution, and is a good representation of group variability in true fMRI datasets. The activation magnitudes were chosen to correspond to a mean SNR of either -10dB or -6dB in each dataset, which corresponds to a realistic SNR in fMRI datasets. In order to model intersubject differences in the SNR, we let this value fluctuate across subjects in the ranges [-18dB -6dB] and [-11 -3dB] respectively. Once again, this kind of fluctuation is a reasonable model of standard inter-subject variability. The CR maps, as well as the RFX maps computed on these datasets, are submitted to Receiver Operating Characteristic (ROC) Analysis. False positive and negative rates were computed while the α parameter (CR) or the threshold (RFX) is varied, resulting in a (sensitivity, specificity) plot. Results were averaged over n = 100 repetitions of the simulation.

  tasks such as passive viewing of horizontal or vertical checkerboards, left or right click after audio or video instruction, computation (subtraction) after video or audio instruction, sentence listening and reading. Events occurred randomly in time (mean inter stimulus interval: 3s), with ten occurrences per event type. 102 right-handed subjects participated in the study. The subjects gave informed consent and the protocol was approved by the local ethics committee. Functional images were acquired on a 3T Bruker scanner using an EPI sequence (TR = 2400ms, TE = 60ms, matrix size=64 × 64, FOV = 24cm × 24cm). Each volume consisted of 34 4mm-thick axial contiguous slices. A session comprised 130 scans. Anatomical T1 images were acquired on the same scanner, with a spatial resolution of 1 × 1 × 1.2 mm 3 . Finally, the cognitive performance of the subjects was controlled using a battery of syntactic and computation tasks. fMRI data pre-processing consisted of 1) temporal Fourier interpolation to correct for between-slice timing, 2) motion estimation; for all subjects, motion estimates were smaller than 1mm and 1 degree, 3) anatomofunctional image coregistration and spatial normalization of the functional images in the MNI/Talairach space.

  here.] Such a result, formulated in terms of regions, readily indicates possible asymmetries in the spatial repartition of activations across subjects: activation in the Inferior Frontal cortex are found in the right hemisphere only, activations in the Supramarginal cortex, the Angular cortex are found in the right hemisphere only and activations in the Precentral regions are more systematic in the right hemisphere (4 regions) than in the left hemisphere (1 region).

Finally, we found

  an almost significant (p < 0.06, after correction for multiple comparisons) ef-fect for the 3D rotation score in the activity in a sub-region the left occipito-parietal boundary [z = 4.03, p < 3.10 -5 , uncorrected at (-33, -65, 45)mm]. Using the ROI-based analysis, the result was also present, and significant [z = 3.10, p < 0.05, corrected at (-27, -73, 36)mm], while this -and only this-ROI had a significant effect of the score on its position across subjects [χ 2 3 = 34.5, p < 0.05, corrected at (-27, -73, 36)mm]. The regions that exhibit significant correlation of their signal or their position with the subject's sex and score in the 3D task are shown in Fig. 11 (a) and (b) respectively.

Fig. 3 .

 3 Fig.3. Illustration of the use of the Belief Propagation algorithm to find correspondences between maxima within a pair of subjects. This is a toy dataset, in a one-dimensional space. The activity maps of subjects s 1 and s 2 are shown on the top of the figure, together with a watershed segmentation. In that case, t s 1 i = 0, t s 1 j = 1 and t s 1 k = 2, while t s 2 l = 0.7, t s 2 m = 1.7 and t s 2 n = 2.7; δτ = 1.4. The related graphs G s 1 and G s 2 are shown below; in this case they are isomorphic. The associations are initialized using Eq. (5), and then refined using Eq. (6): Clearly the message passing algorithm enhances the probabilities P (a s 2 l ← a s 1 i ), P (a s 2 m ← a s 1 j ), and P (a s 2 n ← a s 1 l ), thus compensates the effect of the global translation between the two datasets. This effect is also present with the converse probabilities P (a s 1 . ← a s 2 . ).

Fig. 4 .

 4 Fig.[START_REF] Talairach | Co-Planar Stereotaxic Atlas of the Human Brain. 3-Dimensional Proportional System : An Approach to Cerebral Imaging[END_REF]. Number of false detections obtained with our method in noise only environments (box and whisker plot), compared with its expected value (continuous line), for values of α ranging from 0.1 to 0.8 in 0.05 steps. This is based on 20 runs of 100 simulations under the null hypothesis. The fact that the box plots lie beneath the line shows that the threshold is rather conservative, especially for large values (p > 0.2).

Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 6789 Fig.6. Comparison of the a CR map and RFX maps obtained for a functional contrast that shows regions involved in computation task. This is based on 10 subjects, with one session per subject. (a) The CR map, significant at p < 0.05, corrected level, shows q = 19 active regions; (b) the RFX map is thresholded at the voxel level at p < 10 -3 level, uncorrected ( at p < 0.05 corrected, the map is empty); (c) the RFX map is thresholded at the cluster level at p < 0.05 level, corrected. The CR map clearly shows symmetric parietal regions involved in the computation task, while these regions are not detected with the RFX procedure.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Unsupervised classification of the group of subjects, based on their profile. The dendrogram shows the organization of the population in terms of hierarchical clustering, based on an average linkage approach. The dendrogram shows that there is one main group, plus a few scattered subjects (on the right side).

  Fig. 2. Modeling the density of Ds in the volume. Randomly reshuffling the position of the maxima of activity in subjects σ ∈ {1, .., S} -{s} yields an empirical histogram of Ds under the null hypothesis (green), which can be used to define critical values uα, which can be corrected for multiple comparisons or not, to threshold the density D.

	Subject s 1		Subject s 2
	i	j	k		l	m	n
	0.35
						true histogram
						surrogate histogram
	0.30
	0.25 G s 1		p<0.2, uncorrected	G s 2
	0.20 i	j		k	l	m	n
			p<0.05, corrected
	0.15
	Initialization
	0.05 0.10	P (. ← a s 1 i ) P (. ← a s 1 j ) P (. ← a s 1 k )	P (a s 2 l ← .) P (a s 2 m ← .) P (a s 2 n ← .) 0.59 0.31 0.10 0.42 0.38 0.20 0.25 0.39 0.35
	0 After convergence 0.00	5	10	15	20	25
		P (. ← a s 1 i ) P (. ← a s 1 j ) P (. ← a s 1 k )	P (a s 2 l ← .) P (a s 2 m ← .) P (a s 2 n ← .) 0.74 0.23 0.03 0.41 0.46 0.13 0.14 0.40 0.46

Note that the population is quite homogeneous, mean age=23.9 years and std=3.8 years. In this condition, the absence of an age effect at the group level is not surprising. May

[START_REF] Mériaux | Combined permutation test and mixed-effect model for group average analysis in fMRI[END_REF]

May 29, 2007 DRAFT

MNI coordinates

Anatomical location hemisphere number of subjects (/102) (-54, -55, -10) Temporal Inf R 40 ( 52,[START_REF] Collins | Non-linear cerebral registration with sulcal constraints[END_REF][START_REF] Brett | The problem of functional localization in the human brain[END_REF] Frontal Inf Oper L 51 (-42, 35, 16) Frontal