Alan Tucholka

Virgile Fritsch

Jean-Baptiste Poline

Bertrand Thirion

An empirical comparison of surface-based and volume-based group studies in neuroimaging

Keywords: Cortical surface mapping, fMRI, statistical inference, cluster-level tests, random eects analysis, mixed eects analysis, between-subject variability

come

Introduction

Studying the localization and variability of brain activity across subjects are two crucial aspects of neuroimaging data analysis. This is important both for building models of brain activation organization and to relate inter subjects activity variations to factors of interest, such as behavioral and genetic measurements. In particular, many neuroscience results depend on the the precise localization of the Blood Oxygen-Level Dependent (BOLD) signal changes measured by functional Magnetic Resonance Imaging (fMRI). Controlling the variability of activation position across individuals is essential in order to grant enough sensitivity in activation detection, and to ensure that one-sample group-level inference, that deals with between-subject mean activation, is indeed representative of the true population pattern. Both co-registration and activity detection are most commonly performed in 3D space, by applying linear and non-linear warps to match the anatomical and functional images to a common template, often chosen to be the MNI template [START_REF] Mazziotta | A probabilistic atlas of the humain brain: theory and rationale for its development[END_REF]. Brain activations are then described through the positions of local maxima of activation maps in template space, and these locations are then interpreted in terms of brain regions using standard atlases (see e.g. [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] Toga | Towards multimodal atlases of the human brain[END_REF]).

Activity localization. Once the data are in the common template space, the only information available on the localization of the BOLD signal are the coordinates (x,y,z) in this space. This view precludes a deep understanding of brain organization, which should consider the details of anatomical structure in order to dene properly activity localization [START_REF] Devlin | In praise of tedious anatomy[END_REF]; in general, functional regions should be characterized by their relative position with respect to anatomical landmarks, their extent or their connections to other regions [START_REF] Brett | The problem of functional localization in the human brain[END_REF]. Imperfect spatial registration procedures probably induce an additional blur that is detrimental to the accuracy of brain mapping and to sensitivity.

Activity detection. The standard approach to activation detection [START_REF] Friston | Human Brain Function, chapter 2[END_REF] consists in comparing the images from the dierent subjects on a voxel-by-voxel basis, and to compute a statistical map to test the presence of an activation in each voxel of the standard space. The ensuing multiple testing problem can be addressed directly at the voxel-level [START_REF] Nichols | Controlling the Familywise Error Rate in Functional Neuroimaging: A Comparative Review[END_REF] or by testing the size of clusters dened above a user-chosen threshold [START_REF] Hayasaka | Validating Cluster Size Inference: Random Field and Permutation Methods[END_REF].

Volume-and surface-based spatial normalization

Spatial normalization is therefore crucial to dene accurately the position of functional regions, and thus to detect positive activation at the population level. Many works have been carried out to improve these volume-based co-registration procedures (see [START_REF] Klein | Evaluation of 14 nonlinear deformation algorithms applied to human brain mri registration[END_REF] for a review of most of these procedures). Yet, as a very large part of the data originates from the cortex, methods that work on the cortical surface may be more sensitive than those using the full brain volume. It is well known that volume-based normalization may introduce inaccuracies in anatomical positioning of functional data, estimated up to 1cm in several cortical regions [START_REF] Hellier | Retrospective evaluation of intersubject brain registration[END_REF][START_REF] Stiers | Mapping multiple visual areas in the human brain with a short fMRI sequence[END_REF]. For instance, it is dicult to account for the inter-subject variability of gyri size, shape or position in a 3D referential and such dierences may correspond to the displacement of a functional focus to a dierent gyrus in some subjects.

Surface-based spatial normalization proceeds from a dierent perspective: it consists in the denition of coordinate systems on the cortical surface that match corresponding regions across individuals. This is done by using geometric descriptors of the surface, such as the main sulci [START_REF] Essen | A population-average, landmark-and surface-based (pals) atlas of human cerebral cortex[END_REF][START_REF] Clouchoux | Anatomically constrained surface parameterization for cortical localization[END_REF] or smoothed curvature maps [START_REF] Fischl | Cortical surface-based analysis II: ination, attening and a surface-based coordinate system[END_REF]. A review of these techniques has recently been performed in [START_REF] Klein | Evaluation of volume-based and surface-based brain image registration methods[END_REF]. Surface-based methods do not guarantee that a perfect match will be found across individuals, in particular because sulcal patterns are known to be partly inconsistent even in standard populations [START_REF] Rivière | Automatic recognition of cortical sulci of the human brain using a congregation of neural networks[END_REF]. However, these surface-based methods are probably more accurate for aligning cortical folds than volume-based methods, and they have been shown to correctly align cyto-architectonic boundaries [START_REF] Fischl | Cortical folding patterns and predicting cytoarchitecture[END_REF]. Several studies have shown that a co-registration algorithm based on the cortical surface may better align the functional signal across subjects [START_REF] Fischl | Cortical surface-based analysis II: ination, attening and a surface-based coordinate system[END_REF][START_REF] Argall | Simplied intersubject averaging on the cortical surface using suma[END_REF][START_REF] Desai | Volumetric vs. surface-based alignment for localization of auditory cortex activation[END_REF], 1], but none of these studies used standard group-level random-eects inference procedures, such as a one-sample t-test across individuals.

More generally, besides the fact that non-cortical regions are not analyzed by a surface-based procedure, this analysis suers from an important limitation: it requires that functional data has been accurately projected on the cortical surface, i.e. that the correspondence of functional data with anatomical data is correct. This is a real challenge in practice, as the cortical thickness, but also the distance between the banks of a given sulci or the width of a gyrus are small, in view of the typical resolution of fMRI data [START_REF] Operto | Projection of fmri data onto the cortical surface using anatomically-informed convolution kernels[END_REF]. Moreover, fMRI data are artefacted by geometric distortions related to EPI acquisition. These distortion are not systematically and fully corrected during pre-processing [START_REF] Hutton | Image distortion correction in fmri: A quantitative evaluation[END_REF]. In the present work, we do not focus on the evaluation of the main steps that are necessary to perform adequate surface-based analysis, such as the comparison of functional data projection techniques or EPI distortion correction, but propose dierent methods to compare the results of various group-level inference procedures applied on the cortical surface and in the volume. For these comparisons, we use state-of-the-art group analysis statistical approaches. A procedure has recently been proposed in [START_REF] Greve | Accurate and robust brain image alignment using boundary-based registration[END_REF] to assess the quality of the co-registration between a cortical mesh and fMRI data, and is expected to optimize the spatial match between anatomical structures and functional information. However, to the best of our knowledge, there has been no recent comprehensive investigation on the impact of 2D versus 3D activation detection on fMRI group analysis.

An empirical comparison of volume-and surface-based alignment for functional brain mapping

In this paper, we investigate whether surface-based approaches, that rely on a cortical surface referential, provide better constraints about the position of functional activity, and more precisely, whether this is reected in state-of-the-art inter-subject statistical procedures. Following [START_REF] Fischl | Cortical surface-based analysis II: ination, attening and a surface-based coordinate system[END_REF], we perform functional analysis on the cortical surface for a group of 25 subjects. The inter-subjects analysis relies on matching the subjects cortical surface [START_REF] Fischl | Cortical surface-based analysis II: ination, attening and a surface-based coordinate system[END_REF][START_REF] Fischl | High-resolution intersubject averaging and a coordinate system for the cortical surface[END_REF]. Then we systematically compare the results of surface-and volume-based statistical analysis and provide results on the dierence in sensitivity of the two approaches for dierent tests, for a given control of the type I error. More specically, we use for the comparison mixed-and random-eects inference at the voxel and at the cluster level [START_REF] Mériaux | Combined permutation test and mixed-eect model for group average analysis in fMRI[END_REF] and assess their reproducibility with bootstrap.

This raises a dicult question: how should we compare group results obtained in the volume to those obtained on the cortical surface? In the present work we explore and detail three possibilities, that are based either on the relative activated volume, or the choice of surface-or volume-based spatial referential to compare the results of either kind of analysis. We also address the impact of choosing a surface-or a volume-based smoothing kernel on the results, as this is known to have an important eect on detection and localization in group analysis. We consider an fMRI protocol designed to map quickly several cognitive functions [START_REF] Pinel | Fast reproducible identication and large-scale databasing of individual functional cognitive networks[END_REF] from which we use several contrasts to compare the surface and volume analyses, and provide to cognitive neuroscientists comprehensive information on the sensitivity of those analyses.

The remainder of this paper is organized as follows: in section 2 we rst present the data used in our experiments, and then the processing and inference procedures tested; we describe the results of these tests in section 3 and discuss our results in section 4.

Materials and methods

Dataset used in our experiments

Data were acquired from 25 subjects who performed a functional localizer protocol as described in [START_REF] Pinel | Fast reproducible identication and large-scale databasing of individual functional cognitive networks[END_REF]. 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 dierent 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. Functional images were acquired on an 1.5T General Electric Signa System scanner (General Electric Medical Systems, Milwaukee, WI, USA) using an EPI sequence (TR = 2400ms, TE = 60ms, matrix size = 64×64, FOV = 24cm×24cm, echo spacing = 608µs). Each volume consisted of 40 3mm-thick axial slices without gap. A session comprised 132 EPI scans, of which the rst 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.2 mm, a eld of view of 24 cm and an acquisition matrix of 256 × 256 × 124 voxels, resulting in 124 contiguous double-echo slices with voxel dimensions of (0.9375 × 0.9375 × 1.2) mm 3 .

Pre-processing

The functional data were rst corrected from the EPI distortions using eld maps and the SPM toolbox; the ensuing image resolution was 3.75 × 3.75 × 3mm 3 . Note that the fMRI were considered as initially aligned with the B0 map, and thus were not re-interpolated prior to distortion correction; the application of the correction in itself induced a re-interpolation of the data. Next, a standard pre-processing (correction of dierences in slice timing, rigid-body motion correction and anatomo-functional co-registration) was performed using the SPM5 software on all subjects.

The FreeSurfer software [START_REF] Fischl | Cortical surface-based analysis II: ination, attening and a surface-based coordinate system[END_REF][START_REF] Dale | Cortical surface-based analysis. i. segmentation and surface reconstruction[END_REF], version 4.0.2 was used to segment and reconstruct the cortical surface from T1 MRI data of each subject, and obtained the white matter mesh for both hemispheres. This provides a common spherical coordinate system for each hemisphere in each subject. Pre-processing of the data includes i) segmentation of the white matter, yielding triangular meshes for grey-white and grey-csf (cerebro-spinal uid) interfaces, ii) detection of the deepest sulci, iii) ination of the white surface onto a sphere, iv) deformation to match the deepest sulci positions on the template model.

All data are then converted to the standard GIFTI format for further processing. In order to obtain a node-by-node correspondence, the mesh of each subject was then resampled : i) a regular sphere (icosphere) of diameter equal to the brain-sphere with a reduced number of nodes (about 40k nodes) was created, ii) this sphere was refolded onto the original cortical surface of each subject while preserving node-to-node correspondence with the original icosphere mesh. This corresponds to resampling the individual meshes to match a template mesh, as is usually done in volume-based normalization procedures. Downsampling to 40k nodes entails a minor loss of resolution when compared to fMRI resolution, and provides a four-fold speed gain in the ensuing statistical procedures with respect to the full mesh resolution. An average cortical model across 25 subjects was created for visualization of the results. Each resampled mesh was embedded into the MNI space, then an average brain was obtained by computing the mean 3D position of each node through all subjects in that space. It is important to notice that this mesh is used only for visualization purpose and not in the actual statistical analyses.

Functional images were then projected onto the resampled gray/white interface mesh of each subject using the method described in [START_REF] Operto | Projection of fmri data onto the cortical surface using anatomically-informed convolution kernels[END_REF]. This projection relies on the geometry of the local anatomy: for each node of the mesh of the grey/white interface, a projection kernel is computed, so that it yields a weighted average of the neighboring voxels, where the weights decrease sharply along the normal direction outside of the cortical ribbon, and more smoothly in the tangent plane at the surface. The projection is then performed by convolution of the functional 3D map with the kernels of each node. The parameters of the decay are 5mm in the tangential direction and 2mm in the normal direction. For all subjects and hemispheres, we visually checked that there was no spatial mismatch between the brain mesh and the functional volume. A General Linear Model (GLM) analysis was applied for the volume-and surface-based data using the same analysis code, namely the nipy package http://nipy.sourceforge.net/. The model included the ten conditions of the experiments convolved with a standard hemodynamic lter and its time derivative, a high-pass lter (cuto:128s) and the procedure included an estimation of the noise auto-correlation using an AR(1) model [START_REF] Roche | Incremental activation detection in fMRI series using kalman ltering[END_REF]. Activation maps were derived for the following six functional contrasts (we give short names in italic): i) left-right: left versus right button presses, ii) right-left: right versus left button presses, iii) audio-video: sentence listening versus sentence reading, iv) video-audio: sentence reading versus sentence listening, v) computation-sentences: computation versus sentence reading, vi) reading-visual: reading versus passive checkerboard viewing. These contrasts are expected to reveal dierent aspects of the functional organization of the brain, at dierent localizations, and show dierent sensitivity and reproducibility levels (see [START_REF] Thirion | Analysis of a large fMRI cohort: Statistical and methodological issues for group analyses[END_REF] for a detailed discussion).

In parallel, the fMRI data was also analyzed, using the same linear model, after non-linear spatial normalization, which was performed with the SPM5 software using the new segment method.

This normalization was based on the coregistration with the T1 image that had been warped to the SPM T1 template. Before the linear model application, we used dierent levels of smoothing on the cortical surface and in the volume, corresponding to 0, 4, 8 and 12mm full width at half maximum (fwhm) for the 3D data and 0 or 8mm on the surface, in order to check the impact of smoothing kernel. We also segmented the grey matter in the individual anatomy using SPM5, and created a population-level mask of the voxels that belong to the grey matter with a probability greater than .5; this masked was resampled to the fMRI resolution.

Experiments 2.3.1. Statistical model for group analysis

In this work, S = 25 subjects are considered. For each subject i, and at any location of the domain under consideration (cortical mesh or brain volume), let βi be the estimation of the BOLD eect related to some functional contrast of interest (for notational simplicity, we drop explicit references to the contrast under consideration). βi is distributed around the true eect β i : βi = β i + e i with e i ∼ N (0, s 2 i ) where the estimation variance s 2 i is assumed to be known from the rst-level General Linear Model (GLM) 1 . As is traditionally done in neuroimaging, we further assume that β i = β G + ε i , where ε i ∼ N (0, σ 2 ) and β G is the population-level eect [START_REF] Roche | Mixed-eects statistics for group analysis in fMRI: A nonparametric maximum likelihood approach[END_REF]. We thus have :

βi = β G + ε i , ε i ∼ N (0, σ 2 + s 2 i ), (1)
where σ 2 is the between-subject variance. This is a generalization of the Random Eects (RFX) model in [START_REF] Friston | Human Brain Function, chapter 2[END_REF] which neglects the βi estimation variance, i.e. it assumes s 2 i ≡ 0. Both β G and σ 2 are then estimated by maximizing the log-likelihood of the model specied in Eq. (1) using the Expectation-Maximization (EM) algorithm in [START_REF] Keller | Dealing with spatial normalization errors in fMRI group inference using hierarchical modeling[END_REF]. The following log-likelihood ratios are computed to test the positivity of β G :

L M F X = log sup σ 2 ,β G >0 S i=1 N ( βi ; β G , σ 2 + s 2 i ) sup σ 2 S i=1 N ( βi ; 0, σ 2 + s 2 i ) , (2)
L RF X = log sup σ 2 ,β G >0 S i=1 N ( βi ; β G , σ 2 ) sup σ 2 S i=1 N ( βi ; 0, σ 2 ) (3)
Note that computing L RF X is equivalent to performing a t-test.

Statistical calibration

The distribution of the statistics in Eq. (2-3) under the null hypothesis (β G = 0) is unknown, but can be estimated very simply through a randomization procedure, in which the statistics are recomputed after a sign swap of the observed eects βi . Under the hypothesis that the distribution of the true eects is symmetric about zero everywhere in the brain which is our null hypothesis this procedure yields an exact (possibly conservative) specicity for the test. In order to control the family-wise error rate (FWER), i.e. the probability of detecting one false positive region over the search domain, we consider the distribution of the maximal statistic under the null hypothesis. For a chosen FWER α, this yields a voxel-or vertex-level corrected threshold.

A more sensitive approach to detect extended regions consists in rst thresholding the statistics map at a given level (corresponding e.g. to p < 10 -3 uncorrected), and then to estimate the distribution of size (area (in mm 2 ) or volume (in mm 3 )) of the supra-threshold clusters under the null hypothesis. To solve the multiple comparison issue, the size of the maximal cluster is tabulated under the null hypothesis. Once again, the quantile α of this simulated distribution yields a cluster-level corrected threshold.

In summary, we use the following statistics: voxel-level/vertex-level random eects test (VRFX) and mixed eects test (VMFX), cluster-level random eects test (CRFX) and mixed eects test (CMFX). We will always refer to surface-based or volume-based statistics to make the underlying spatial model explicit.

Comparison of surface-based and volume-based results

Next we discuss how to compare the dierence in sensitivity between the spatial-and volumebased analyses. One rst challenge is that this comparison is performed in the absence of a ground truth, so that it is dicult to characterize the correctness of surface-or volume-domain detections. However, the procedure that we use has the theoretical guarantee of providing the same rate of false detections (type I error) for each test. This holds separately for univariate tests on the one hand, and cluster size tests on the other hand (the sensitivity of both tests are not easily compared directly as they do not have the same regional specicity). In these conditions, the fact that a spatial procedure yields the detection of more regions means that this procedure is better suited to detect the functional activity related to this contrast. The diculty is thus to quantify the relative amount of detections performed in either domain, as the results are observed and characterized in two dierent spaces. In this work, we describe and perform three possible comparisons:

• The rst one is to consider the test sensitivity, i.e. the relative amount of activated regions in either domain (surface or volume). This can give an indication on the performance of either domain to outline active regions, but this is confounded by the fact that 3D maps include sub-cortical regions, in which activations may (cerebellum, deep nuclei) or may not (white matter) be found. We have tried to remove this confound by running also the volume-based analysis on a mask of the grey matter.

• A second possible comparison is to project the result of group-level volume-based analyses onto the cortical surface of each subject, and then to report the activated surface (in mm 2 ) obtained, as compared to the surface declared active after full surface-based analysis. In the absence of a satisfactory population-level model of the cortical surface, we iteratively took each individual surface as reference space, and report the mean and standard deviation of the activated surface across these models. This avoids the confounds previously described, but it is dicult to decide which vertices are active when we recombine the projection of the group results on a given individual anatomy: in this work, we decided to declare active all nodes that reached a signal level of .5 after projection of a binary activation maps (0 for inactive, 1 for active voxels) onto the cortical surface, corresponding to a 50% condence that the node is indeed active.

• Finally, a third solution consists in embedding the activated regions obtained on the cortical surface into the MNI space, and then to compare the two sets of positions in MNI space. The position comparison procedure is dened using a kernel-base metric that measures the discrepancy between dierent sets of positions, and is described in Section 2.4.1. Again, as the embedding depends on the chosen cortical surface, we used the grey/white interface mesh of each subject and average the result across all 25 subjects. The embedding and projection approaches are illustrated in Fig. 1. Finally, we also considered the impact of choosing the medial cortical surface instead of the grey-white matter interface in the projection approach.

Kernel similarity between active regions

Here we dene the similarity measure used to compare the coordinates obtained from surfaceand volume-based analysis in the MNI space: Let t surf = (t surf i ) i=1..N and t vol = (t vol i ) i=1..V be the set of coordinates of supra-threshold positions obtained from surface-based and volume-based analysis respectively, N and V being the number of nodes and voxels under consideration. We dene the two quantities:

ψ(t vol ; t surf ) = mean v=1..V max n=1..N exp - t vol v -t surf n 2 2δ 2 (4)
Embedding approach (a)

Projection approach and

(d) (d) (b) (c) (b) Subject s ∈[1..S] Subject s ∈[1..S]
ψ(t surf ; t vol ) = mean n=1..N max v=1..V exp - t vol v -t surf n 2 2δ 2 , ( 5
)
where δ is a xed distance (δ = 6mm in our experiments Regarding the choice of δ, the meaning of this parameter is that, whenever the distance between one node and it nearest voxel is larger that 2 × δ, the contribution to ψ is close to 0, meaning that the node has no corresponding active voxel in the volume. δ measures the admissible discrepancy of the position of activations in MNI space, and we found δ = 6mm to be a reasonable choice. However, our evaluation metric varies smoothly with respect to this parameter, so that the results are stable with respect to this choice too.

Bootstrap reproducibility of the active regions

Finally, we performed an analysis of the reproducibility of the active regions in 2D and 3D using a bootstrap procedure: we draw S subject with replacement from the initial population B times, perform the statistical analysis, and compute a replication map R B which yields the counts of activation of each site. A reproducibility index is derived below from this replication map. Although a measure based on a binomial mixture has been proposed previously for this purpose [START_REF] Genovese | Estimating test-retest reliability in functional MR imaging. I: Statistical methodology[END_REF][START_REF] Maitra | Test-retest reliability estimation of functional MRI data[END_REF][START_REF] Thirion | Analysis of a large fMRI cohort: Statistical and methodological issues for group analyses[END_REF], we use here a direct estimation of the conditional probability ρ that a supra threshold voxel in one particular dataset will also be above threshold in another sample. This was simply obtained by counting the average number of co-activations across pairs of experiments, divided by the average number of activation detections across experiments. Let h be the histogram of R B : (h b ) 1≤b≤B counts the number of times a voxel/node has been declared active among B replications; then let

ρ = 1 B -1 B b=1 b(b -1)h b B b=1 bh b (6)
ρ ∈ [0, 1] measures the probability that a given active site will in average, be reproduced in a following experiment. In our experiments, we use B = 10 and S = 25.

For visualization, we used the Anatomist software.

Results

A summary of the active regions activations found for the six contrasts under investigation are given in Fig. 2.

We successively report the results of the following tests for six dierent functional contrasts studied in the protocol: • Comparison of the relative amount of activated regions found in surface-and volume-based analysis.

• Comparison of surface-based and volume-based supra-threshold regions through kernel-based similarity.

• Reproducibility of the supra-threshold regions in the surface of in the volume.

The level of smoothing can have a confounding eect on the conclusions that might be drawn from these experiments. We thus systematically considered dierent levels of smoothing: 0, 4, 8 and 12 mm fwhm in the volume, 0 or 8 mm on the cortical surface. We report all or part of these results in the following experiments. Surface-based smoothing was performed on the mesh, using a tool from the Brainvisa toolbox http://brainvisa.info/index.html.

Relative sensitivity 3.1.1. Rate of signicantly active regions

The proportion of supra-threshold surface or volume is given in Table 1 for mixed-eects inference, and in supplementary material for random-eects inference (see Table A.3). We observe that: i) cluster-level statistics systematically yield more extended supra-threshold regions than voxel-or node-level statistics; ii) mixed-eects statistics are more sensitive than random-eects statistics; iii) surface-based analysis systematically yields a larger relative portion of activated regions than volume-based analysis. Results i) and ii) are common ndings, while iii) provides a rst empirical conrmation of the advantage of surface-based analyses.

In the sequel, we no longer report results with VRFX/CRFX, since they are qualitatively similar to those obtained with VMFX and CMFX respectively, with a general tendency toward less sensitivity and reproducibility.

Next, we consider the impact of the level of smoothing on these quantities; to simplify the presentation, we only consider here a smoothing kernel of 8mm fwhm, applied on the surface or in the volume. The results are given in Table 2 using mixed-eects inference at the voxel/node (VMFX) level only. We also include the result in the volume, but limited to the grey matter, without smoothing. This should be compared with the results provided in Table 1, that gives corresponding results without smoothing.

We observe that the tendency observed previously still holds after smoothing the data: the proportion of activations is larger on the surface, and that smoothing increases this proportion. Restricting the volume analysis to grey matter increases the values, but it still remains below surface-based results.

Spatial extent of signicantly active regions on the cortical surface (projection approach)

The spatial extent of the supra-threshold regions after surface-based analysis is compared with the extent of the 3D active regions of the volume-based group map, reprojected onto each individual surface. The average activated region is reported, as explained in section 2.4. The results are shown in Fig. 3 for the six contrasts under study, using the voxel-level mixed eects statistics (VMFX). We also consider the impact of smoothing the data in the volume or on the surface. We also ran the experiment on a mask of the grey matter in the volume, but did not nd any dierence with respect to a full-volume study.

As expected, we clearly see that surface-based analysis yields a larger portion of the cortical surface, and that smoothing the data tends to increase the relative size of active regions. It is also important to note that the amount of activated surface after re-projection depends only weakly on the individual anatomy used for projection, as illustrated by the relatively low variability of the area of the projection.

Kernel-based sensitivity analysis (embedding approach)

Next, we investigate the similarity of the surface and volume-based patterns by considering the natural embedding the cortical surface in MNI space, and then assessing the similarity of the positions of active regions using the metric detailed in section 2.4.1. The results are shown in Fig. 4 for six contrasts, using CMFX and VMFX inference, and with two levels of smoothing: 0mm and 8mm.

We observe a very high similarity overall, with similarity values close to .8, meaning that on average, an active voxel can be found within d = 4mm ((expd 2 2δ 2 = .8) yields d 4mm, given that δ = 6mm) of an active node, and conversely. However, some dierences can be noted, and in particular between voxel-based and cluster-based inference: without additional smoothing ψ(t surf ; t vol ) is similar to ψ(t vol ; t surf ) in cluster-level inference, while ψ(t vol ; t surf ) > ψ(t surf ; t vol ) in voxel/vertex-level inference; when an 8mm smoothing kernel is applied to volume and surface data, ψ(t surf ; t vol ) > ψ(t vol ; t surf ) in cluster-level inference, while ψ(t vol ; t surf ) ψ(t surf ; t vol ) in voxel/vertex-level inference. This readily means that, in peak-level inference, the surface-based pattern seems to reveal more details (activation peaks) than voxel-based analysis. This eect vanishes in cluster-level inference, or when the level of smoothing is increased. The combination of large smoothing and cluster-level inference shows regions in the volume that are not matched with surface-based analysis. Note that we also investigated the use of the medial cortical surface of the grey-white matter interface when embedding the surface activation into individual space, but the dierence with respect to the use of the grey-white matter interface is barely noticeable, as the distance between the two surfaces is much smaller than δ. We also considered using the activations found in greymatter volume only, but the eect is also negligible.

Qualitative comparison

In order to compare more qualitatively the output of surface-based versus volume-based inference, we present both types of results on the same gure in Fig. 5. We present results for mixed eects analysis, given that random-eects yield very similar results, albeit with less sensitivity. We present volume-based and surface-based inference results for two contrasts, computation-sentences and reading-visual.

With the computation-sentences contrast, some active surface portions are relatively far from active voxels (but not the converse) in voxel-level inference (Fig. 5, rows 1 and 3), but not in cluster-level inference (Fig. 5, rows 2 and 4). This is consistent with the fact that ψ(t vol ; t surf ) > ψ(t surf ; t vol ) in peak-level inference. This is one of the cases, where surface-based analysis can be viewed as more sensitive than volume-based analysis. This eect is less obvious with the reading-visual contrast.

In cluster level inference, we obtain a greater consistency, i.e. all the supra-threshold locations in either domain correspond to supra-threshold location in the other domain. The most striking result is that, in the case of reading-visual contrast, volume-based analysis yields a huge cluster that gathers the superior temporal, Broca's and pre-central gyri, while these appear as distinct clusters in surface-based analysis. This eect can be considered as problematic, as the huge volume-based region clearly results from the agglomeration of spatially distant active foci observed in dierent subjects. Our interpretation is that the spatial variability of these foci creates a very large region where the average signal across subjects is positive; when many subjects are observed (here S = 25), the mean deviation to zero becomes statistically signicant, and a very large supra-threshold cluster appears. Surface-based analysis nicely avoids this eect, as it produced spatially distant active regions on the cortical surface.

Reproducibility analysis

Finally the reproducibility statistic has been computed, based on a bootstrap procedure described in section 2.5. The results are given in Fig. 6 in the case of cluster-level or voxel-level inference, using the MFX statistic, for the six contrasts of interest.

First of all, the reproducibility values are in the [. [START_REF] Clouchoux | Anatomically constrained surface parameterization for cortical localization[END_REF], 1] range, reecting that there is a relatively high probability of reproducing the results at the voxel/vertex level when sampling a new group in the same population; this overall reproducibility level depends on the contrast under investigation, as motor and auditory activations yield the most reproducible activity patterns. As could be anticipated, smoothing increases the level of reproducibility. It can be seen that in all cases, surface-based analysis yields higher reproducibility scores than volume-based inference. There is a simple geometric reason for this: projecting the data from 3D to 2D removes one dimension of data variability (the dimension that is perpendicular to the surface), so that the apparent and the reading-visual contrasts (rows 3 and 4), using voxel-level (rows 1 and 3) or cluster-level (rows 2 and 4) inference. Surface activations are plotted as a red-yellow texture on the mesh and active voxels are printed in blue. Note the overall correspondence of both detection procedures in each contrast, and the dissociation between peak statistics, where the surface representation yields a wider network, and cluster size inference, where the volume representation yields extended, but non-specic networks. Maps are thresholded at the p < 0.05, corrected level. cross-subject variance is reduced after projection. Additionally, the better correspondence of cortical structures in the surface-based coordinate system should further decrease the variability of supra-threshold patterns.

Discussion

The results obtained can be summarized as follows. In all our experiments, we observed that:

• Cluster-level statistics systematically yield more extended supra-threshold regions than voxelor node-level statistics. This was expected, as cluster-level statistics that test the size of supra-threshold regions take advantage of the intrinsic smoothness of the data that is ignored by peak statistics [START_REF] Hayasaka | Validating Cluster Size Inference: Random Field and Permutation Methods[END_REF]. This advantage of cluster level inference is mitigated by two important drawbacks: i) it relies on an arbitrary cluster-forming threshold for which there is no straightforward or simple optimal choice (see [START_REF] Smith | Threshold-free cluster enhancement: addressing problems of smoothing, threshold dependence and localisation in cluster inference[END_REF] for a more complete discussion on this topic); ii) it provides a weak control on false detections, i.e. no guarantee on the true status (active or inactive) of each voxel or vertex within supra-threshold regions.

• Mixed-eects statistics are more sensitive than random-eects statistics. Although this increased sensitivity is not guaranteed in theory, it is often observed that mixed eects statistics, that take into account both intra-subject and one inter-subject variance terms, use more information from the input data, and thus allow more sensitive inference [START_REF] Roche | Mixed-eects statistics for group analysis in fMRI: A nonparametric maximum likelihood approach[END_REF]. In particular, these approaches down-weight the observations that are less reliable (i.e. with high rst-level variance), and are thus less sensitive to artefactual values.

• Surface-based analysis systematically yields a larger portion of supra-threshold regions than volume-based analysis. This eect was expected, as surface based analysis focuses by definition on the cortical grey matter, and thus discards many regions where BOLD activity is not expected (but also some regions where BOLD activity might occur). But, as we still obtain a signicant dierence in favor of surface based analysis when restricting the volume to the grey matter, the dierence should rather be attributed to better between-subject correspondences in the surface space.

• An embedding of the cortical surface into each subject's space reveals a dissociation: peaklevel inference shows more details on the surface than in the volume, while cluster-size inference displays larger regions in the volume. This dissociation can be explained as follows:

in peak-level inference, surface based may be more sensitive because the search domain does not include white matter and should yield better inter-subject correspondence. Reciprocally, cluster-level inference in the volume shows very large clusters above a given threshold, because limited spatial accuracy in volume normalization tends to merge a network of regions into a large supra-threshold region: this can be viewed as a spatial jitter eect. This is clearly seen in Fig 5.

• Supra-threshold regions are more stable on the surface than in the volume. This represents another advantage of surface-based inference: the geometric variability of active regions is simply shrunk by dimension reduction, but also by better registration between subjects.

To summarize, the overall greater sensitivity of surface-based analysis is explained by two factors:

• The search domain is smaller (i.e. it consists in the cortex only), while the volume-based approach tests also white matter and cerebro-spinal uid regions. The price to pay is that sub-cortical structures (thalamus, basal ganglia, cerebellum etc.) are not considered in the surface-based representations. However, our experiments show that this does not explain all the observed dierences.

• The co-registration of the data is certainly more accurate, as often discussed in the literature [START_REF] Fischl | Cortical surface-based analysis II: ination, attening and a surface-based coordinate system[END_REF][START_REF] Fischl | Cortical folding patterns and predicting cytoarchitecture[END_REF], but hard to prove in practical settings. Note that this better co-registration could have the eect to reduce the spread of signicant regions by removing spatial uncertainty on their position.

When considering cluster-level statistics, it can be expected that the ner co-registration related to the surface-based approach has a mitigated eect on the detection of large clusters: volumebased approaches tend to join smaller clusters which are close in Euclidean space, but relatively distant on the geodesic (surface-based) representation. It should also be reminded that clusterlevel analyses allow only a weak control of false detections (see e.g. [START_REF] Ge | Resampling-based multiple testing for microarray data analysis[END_REF]): a supra-threshold cluster should contain at least one voxel for which the null hypothesis can be rejected, but it cannot be concluded that all the voxels or nodes in that cluster are indeed active.

However, in surface-based approaches, if the resolution is not ne enough, as this is the case in our experiments, or if EPI distortions are not perfectly corrected, the functional signal may be diluted to neighboring gyri, creating false supplementary clusters (as we can see on Fig. 7). The eect might be consistent across subjects and thus results in a spurious activation focus in a group analysis. Note that volume-based analyses suer from the same inaccuracy, but that it appears only when activations are displayed on the cortical surface.

It is important to notice that we did not do any systematic study of dierent pipelines, and that part of the results described here would be slightly altered when using dierent pipelines and pre-processing tools, e.g. using only FreeSurfer pipeline, and that constant improvement of these tools also continuously changes the picture. In particular, the three following points should be considered carefully: i) the number of data resamplings performed along the analysis pipeline (three in the present case: the rst one during distortion correction, a second one during motion correction and the third one during the projection onto the surface), ii) the potential mismatch remaining between EPI and T1 data, in particular when the acquisitions do not entirely cover the brain, e.g. if the tip of the motor cortex is not within the acquired volume, as it was the case in some subjects in our dataset and iii) the choice and parametrization of the projection method, which has to trade o sensitivity, e.g. integrating over a larger domain to avoid missing some signal, against specicity, avoiding bleeding across regions or tissues. We defer a more exhaustive study of these dierent processing steps to future work. Still, we believe that the main eects described here would carry on to many surface-based studies. It is hard to predict whether stronger fMRI contrasts would alter the conclusions: on the one hand, the decrease of false negatives may benet surface-based analyses more than volume-based analyses; on the other hand, the sensitivity may saturate. An interesting question that can still be addressed is how a surface co-registration procedure based on geometric features such as FreeSurfer would compare to representations of the surface that take The activity on the pre-central gyrus for a computation task has also been projected on the neighboring gyrus (post-central gyrus) in a fraction of the subjects, and the ensuing eect can be detected in a group study (p < 0.05 corrected, 25 subjects, VMFX statistic). It is presented on the average cortical surface (left) and on an inated representation (right). This kind of eect can be due to poor initial fMRI data resolution, to the limitations of EPI distortion correction procedures, and to the smoothing induced by successive interpolation steps on the BOLD data. This map was obtained using the projection method in [START_REF] Operto | Projection of fmri data onto the cortical surface using anatomically-informed convolution kernels[END_REF]. (bottom) The limited resolution of standard fMRI data ∼ (3mm) 3 makes this kind of eect hardly avoidable. In this schematic case, the volume-based approach would detect only one large cluster (in blue), while a surface-based approach would project the activated region on both the 2 gyri and would detect 2 dierent clusters. into account sulcus labelling. Some experiments have indeed shown that sulcus-based coordinate systems tend to stabilize the position of some functional landmarks [START_REF] Tucholka | Triangulating cortical functional networks with anatomical landmarks[END_REF]. The combination of shape information exacted from gyral curves, cortical surfaces together with intensity images [START_REF] Postelnicu | Combined volumetric and surface registration[END_REF][START_REF] Du | Whole brain dieomorphic metric mapping via integration of sulcal and gyral curves, cortical surfaces, and images[END_REF] is a promising approach in that respect. Some algorithms also include include information from diusion MRI, such as the position of the main ber bundles [START_REF] Siless | Joint T1 and Brain Fiber Dieomorphic Registration Using the Demons[END_REF].

Conclusion.

Performing fMRI group analysis on the cortical surface instead of the brain volume may benet the detection of some foci of activity, and thus yield more specic, and possibly more sensitive analysis than standard volume-based approaches. It may reveal sharper contrasts in the functional data and thus provide more reliable markers of brain functional anatomy. A cortical surface analysis is quite costly in term of data processing, but it detects at least as many cortical regions as volume-based methods, with better reproducibility, and the additional guarantee of providing more spatially specic clusters. While it seems benecial to group analysis, it is limited by the issues related to potentially incorrect projection of fMRI signal onto the cortical surface, which appears to be the most risky step for surface-based inference. Future work may therefore be directed at providing quality check measures of the accuracy of the projection.

Figure 1 :

 1 Figure1: Illustration of the two main methodologies for the comparison of surface-based and volume-based group analysis results. In the projection approach (top), the activation found in MNI space (d) are projected onto the the cortical surface of each individual (c), and thus compared to the active areas found in the surface-based analysis, then represented in the cortical geometry of the same subject (b). In the embedding approach (bottom) the activations obtained on the surface in the common space represented by a sphere (a), are embedded in the 3D space related to the cortical geometry of each individual in the MNI space (b), and then compared to the volume activation found in MNI space (d).

Figure 2 :

 2 Figure 2: Activations found in the six contrasts under investigation. The results are represented for a mixed-eects model tested with cluster-level inference. The left plots show the active regions on the cortex, for the left and right hemisphere, using inner and outer views. The rights plot shows the volume maps centered at the map-level peak, and all the maps are thresholded at the p < .05 corrected, level.

Figure 3 :

 3 Figure3: Comparison of the extent (in mm 2 ) of supra-threshold regions on the surface after full analysis on the surface or after re-projection of the 3D results onto the surface. The results are shown using the VMFX statistics on the six contrasts and for various smoothing kernel sizes, on the surface and in the volume. Similar trends are observed using other statistics, cluster-level inference, and stronger smoothing. Error bars represent the variability associated with the choice of the individual cortical surface as projection space.

Figure 4 :

 4 Figure 4: Kernel-based similarity of surface-based and volume-based activated regions. The bars in blue indicate how close the position of supra-threshold voxels is to nearby supra-threshold cortical nodes in average; reciprocally, the bars in red indicate how close supra-threshold nodes are to nearby supra-threshold voxels on average. The results are shown for six contrasts, using voxel-level analysis (a, b) and cluster-level (c, d) inference with the mixedeect model, without smoothing (a, c) or after 8mm smoothing (b, d). Dierences at the p < 0.05 uncorrected level, are highlighted with one star, and results signicant at p < 0.001 level are highlighted by two stars. The error bars correspond to the choice of the individual surface chosen to embed active regions into MNI space.

Figure 5 :

 5 Figure 5: Volume-based versus surface-based results comparisons for the computation-sentences (rows 1 and 2)and the reading-visual contrasts (rows 3 and 4), using voxel-level (rows 1 and 3) or cluster-level (rows 2 and 4) inference. Surface activations are plotted as a red-yellow texture on the mesh and active voxels are printed in blue. Note the overall correspondence of both detection procedures in each contrast, and the dissociation between peak statistics, where the surface representation yields a wider network, and cluster size inference, where the volume representation yields extended, but non-specic networks. Maps are thresholded at the p < 0.05, corrected level.

Figure 6 :

 6 Figure6: Reproducibility index for a voxel-/node-level inference on the surface or in the volume, after dierent amounts of smoothing, based on a bootstrap procedure. We can observe that surface-based inference systematically yields more reproducible results than volume-based inference for all levels of smoothing considered.

Figure 7 :

 7 Figure 7: Potential issues with fMRI data projection to the cortical surface. (top)The activity on the pre-central gyrus for a computation task has also been projected on the neighboring gyrus (post-central gyrus) in a fraction of the subjects, and the ensuing eect can be detected in a group study (p < 0.05 corrected, 25 subjects, VMFX statistic). It is presented on the average cortical surface (left) and on an inated representation (right). This kind of eect can be due to poor initial fMRI data resolution, to the limitations of EPI distortion correction procedures, and to the smoothing induced by successive interpolation steps on the BOLD data. This map was obtained using the projection method in[START_REF] Operto | Projection of fmri data onto the cortical surface using anatomically-informed convolution kernels[END_REF]. (bottom) The limited resolution of standard fMRI data ∼ (3mm) 3 makes this kind of eect hardly avoidable. In this schematic case, the volume-based approach would detect only one large cluster (in blue), while a surface-based approach would project the activated region on both the 2 gyri and would detect 2 dierent clusters.

  ). ψ(t surf ; t vol ) and ψ(t vol ; t surf ) quantify how well the values in (t vol ) and (t surf ) correspond to each other: A high value of ψ(t surf ; t vol ) indicates that, for each supra-threshold region on the surface, a close active region can in general be found in the volume. Reciprocally, a high value of ψ(t vol ; t surf ) indicates that each supra-threshold voxel in the volume is closely matched by a supra-threshold node on the surface. Put dierently 1 2 ψ(t vol ; t surf ) + ψ(t surf ; t vol ) ∈ [0, 1] provides a concordance index of the supra-threshold region on both domains, while a positive dierence ψ(t vol ; t surf ) -ψ(t surf ; t vol ) indicates that volumebased activations are more frequently found in the vicinity of surface-based activations than the converse, i.e. that surface-based analysis is more sensitive.

Table 1 :

 1 Relative volume (in percent) of activated areas for dierent statistical inference procedures performed either on the surface or in the volume. The tests are performed for 6 dierent contrasts (in column), and the results are provided with voxel/vertex-level mixed eects inference (VMFX) and cluster-level mixed eects inference (CMFX). The threshold p-value is 0.05, corrected.

	Contrasts	computation-	reading-	left-	right-	audio-	video-
		sentences	visual	right	left	video	audio
	VMFX, volume 0.126	0.975	0.287	0.317	1.018	0.755
	VMFX, surface 1.298	3.571	1.754	1.114	3.562	5.474
	CMFX, volume 2.075	4.172	1.337	1.25	3.193	4.906
	CMFX, surface 9.221	8.96	3.731	2.326	5.907	11.965

Table 2 :

 2 Relative volume (in percent) of activated areas using voxel/vertex-level mixed eects inference (VMFX), on the surface and in the volume, for a smoothing kernel of 8mm fwhm. The tests are performed for 6 dierent contrasts (in column). The last row yields results on volume data, restricted to the grey matter and without smoothing.

	Contrasts	computation-	reading-	left-	right-	audio-	video-
		sentences	visual	right	left	video	audio
	volume, 8mm smoothing	1.312	3.887	1.133	1.287	3.409	3.846
	surface, 8mm smoothing	5.802	6.966	2.9	1.951	5.687	9.602
	GM volume, no smoothing 0.204	1.516	0.447	0.483	1.634	1.166

As there are often many degrees of freedom at the rst level the error made here is likely to be small.

Acknowledgments. The authors acknowledge support from the ANR grant ViMAGINE ANR-08-BLAN-0250-02, BrainPedia ANR-10-JCJC-1408-01 and from the Digiteo DIM grant HiDiNim. We would also like to thank Philippe Pinel for providing the data we used in this experiment.

[1] Anticevic, A., Dierker, D. L., Gillespie, S. K., Repovs, G., Csernansky, J. G., Essen, D. C. V.,