Alexandre Abraham
email: abraham.alexandre@gmail.com

Michael Milham

Adriana Di Martino

R Cameron Craddock

Dimitris Samaras

Bertrand Thirion

Gaël Varoquaux

Adriana Di

Gael Varoquaux

Deriving reproducible biomarkers from multi-site resting-state data: An Autism-based example

Keywords: data heterogeneity, resting-state fMRI, data pipelines, biomarkers, connectome, autism spectrum disorders

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

Introduction

In psychiatry, as in other fields of medicine, both the standardized observation of signs, as well as the symptom profile are critical for diagnosis. However, compared to other fields of medicine, psychiatry lacks accompanying objective markers that could lead to more refined diagnoses and targeted treatment [START_REF] Kapur | Why has it taken so long for biological psychiatry to develop clinical tests and what to do about it&quest[END_REF]. Advances in non-invasive brain imaging techniques and analyses (e. g. [START_REF] Craddock | Imaging human connectomes at the macroscale[END_REF][START_REF] Van Essen | The future of the human connectome[END_REF]) are showing great promise for uncovering patterns of brain structure and function that can be used as objective measures of mental illness. Such neurophenotypes are important for clinical applications such as disease staging, determination of risk prognosis, prediction and monitoring of treatment response, and aid towards diagnosis (e. g. [START_REF] Castellanos | Clinical applications of the functional connectome[END_REF]).

Among the many imaging techniques available, restingstate fMRI (R-fMRI) is a promising candidate to define functional neurophenotypes [START_REF] Kelly | Competition between functional brain networks mediates behavioral variability[END_REF][START_REF] Van Essen | The future of the human connectome[END_REF]. In particular, it is noninvasive and, unlike conventional task-based fMRI, it does not require a constrained experimental setup nor the active and focused participation of the subject. It has been proven to capture interactions between brain regions that may lead to neuropathology diagnostic biomarkers [START_REF] Greicius | Resting-state functional connectivity in neuropsychiatric disorders[END_REF]. Numerous studies have linked variations in brain functional architecture measured from R-fMRI to behavioral traits and mental health conditions such as Alzheimer disease (e. g. [START_REF] Greicius | Default-mode network activity distinguishes Alzheimer's disease from healthy aging: evidence from functional MRI[END_REF][START_REF] Chen | Classification of alzheimer disease, mild cognitive impairment, and normal cognitive status with large-scale network analysis based on resting-state functional mr imaging[END_REF], Schizophrenia (e. g. [START_REF] Garrity | Aberrant" default mode" functional connectivity in schizophrenia[END_REF][START_REF] Zhou | Functional disintegration in paranoid schizophrenia using resting-state fMRI[END_REF][START_REF] Jafri | A method for functional network connectivity among spatially independent resting-state components in schizophrenia[END_REF][START_REF] Calhoun | Exploring the psychosis functional connectome: aberrant intrinsic networks in schizophrenia and bipolar disorder[END_REF]), ADHD, autism (e. g. [START_REF] Plitt | Functional connectivity classification of autism identifies highly predictive brain features but falls short of biomarker standards[END_REF]) and others (e. g. [START_REF] Anderson | Abnormal brain synchrony in down syndrome[END_REF]). Extending these findings, predictive modeling approaches have revealed patterns of brain functional connectivity that could serve as biomarkers for classifying depression (e. g. [START_REF] Craddock | Disease state prediction from resting state functional connectivity[END_REF]), ADHD (e. g. [START_REF] Consortium | The adhd-200 consortium: a model to advance the translational potential of neuroimaging in clinical neuroscience[END_REF]), autism (e. g. [START_REF] Anderson | Functional connectivity magnetic resonance imaging classification of autism[END_REF]), and even age [START_REF] Dosenbach | Prediction of individual brain maturity using fMRI[END_REF]. This growing number of studies has shown the feasibility of using R-fMRI to identify biomarkers. However questions about the readiness of R-fMRI to detect clinically useful biomarkers remain [START_REF] Plitt | Functional connectivity classification of autism identifies highly predictive brain features but falls short of biomarker standards[END_REF]. In particular, the reproducibility and generalizability of these approaches in research or clinical settings are debatable. Given the modest sample size of most R-fMRI studies, the effect of cross-study differences in data acquisition, image processing, and sampling strategies [START_REF] Desmond | Estimating sample size in functional MRI (fMRI) neuroimaging studies: statistical power analyses[END_REF][START_REF] Murphy | An empirical investigation into the number of subjects required for an event-related fmri study[END_REF][START_REF] Thirion | Analysis of a large fMRI cohort: Statistical and methodological issues for group analyses[END_REF] has not been quantified.

Using larger datasets is commonly cited as a solution to challenges in reproducibility and statistical power [START_REF] Button | Power failure: why small sample size undermines the reliability of neuroscience[END_REF]. They are considered a prerequisite to R-fMRIbased classifiers for the detection of psychiatric illness. Recent efforts have accelerated the generation of large databases through sharing and aggregating independent data samples [START_REF] Fair | Distinct neural signatures detected for ADHD subtypes after controlling for micro-movements in resting state functional connectivity MRI data[END_REF][START_REF] Mennes | Making data sharing work: The fcp/indi experience[END_REF][START_REF] Martino | The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism[END_REF]. However, a number of concerns must be addressed before accepting the utility of this approach. Most notably, the many potential sources of uncontrolled variation that can exist across studies and sites, which range from MRI acquisition protocols (e. g. scanner type, imaging sequence, see [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]), to participant instructions (e. g. eyes open vs. closed, see [START_REF] Yan | Standardizing the intrinsic brain: towards robust measurement of inter-individual variation in 1000 functional connectomes[END_REF]), to recruitment strategies (age-group, IQ-range, level of impairment, treatment history and acceptable comorbidities). Such variation in aggregate samples is often viewed as dissuasive, as its effect on diagnosis and biomarker extraction is unknown. It commonly motivates researchers to limit the number of sites included in their analyses at the cost of sample size.

Cross-validated results obtained from predictive models are more robust to inhomogeneities: they measure model generalizability by applying it to unseen data, i. e., data not used to train the model. In particular, leave-out cross-validation strategies, which remove single individuals (or random subsets), are common in biomarkers studies. However, these strategies do not measure the effect of potential site-specific confounds. In the present study we leverage aggregated R-fMRI samples to address this problem. Instead of leaving out random subsamples as test sets, we left out entire sites to measure performance in the presence of uncontrolled variability.

Beyond challenges due to inter-site data heterogeneity, choices in the functional-connectivity data-processing pipeline further add to the variability of results [START_REF] Carp | The secret lives of experiments: methods reporting in the fMRI literature[END_REF][START_REF] Yan | Standardizing the intrinsic brain: towards robust measurement of inter-individual variation in 1000 functional connectomes[END_REF][START_REF] Shirer | Optimization of rs-fMRI pre-processing for enhanced signalnoise separation, test-retest reliability, and group discrimination[END_REF]. While preprocessing procedures are now standard, the different steps of the prediction pipeline vary from one study to another. These entail specifying regions of interest, extracting regional time courses, computing connectivity between regions, and identifying connections that relate to subject's phenotypes [START_REF] Craddock | Disease state prediction from resting state functional connectivity[END_REF][START_REF] Richiardi | Decoding brain states from fMRI connectivity graphs[END_REF][START_REF] Shirer | Decoding subject-driven cognitive states with whole-brain connectivity patterns[END_REF][START_REF] Eickhoff | Connectivity-based parcellation: Critique and implications[END_REF].

Lack of ground truth for brain functional architecture undermines the validation of R-fMRI data-processing pipelines. The use of functional connectivity for individual prediction suggests a natural figure of merit: prediction accuracy. We contribute quantitative evaluations, to help settling down on a parameter-free pipeline for R-fMRI. Using efficient implementations, we were able to evaluate many pipeline options and select the best method to estimate atlases, extract connectivity matrices, and predict phenotypes.

To demonstrate that pipelines to extract R-fMRI neuro-phenotypes can reliably learn inter-site biomarkers of psychiatric status on inhomogeneous data, we analyzed R-fMRI in the Autism Brain Imaging Data Exchange (ABIDE) [START_REF] Martino | The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism[END_REF]. It compiles a dataset of 1112 R-fMRI participants by gathering data from 17 different sites. After preprocessing, we selected 871 to meet quality criteria for MRI and phenotypic information. Our inter-site prediction methodology reproduced conditions found under most clinical settings, by leaving out whole sites and using them as newly seen test sets. To validate the robustness of our approach, we performed nested cross-validation and varied samples per inclusion criteria (e. g. sex, age). Finally, to assess losses in predictive power associated with using a heterogeneous aggregate dataset instead of uniformly defined samples, we included a comparison of intraand inter-site prediction strategies.

Material and methods

A connectome is a functional-connectivity matrix between a set brain regions of interest (ROIs). We call such a set of ROIs an atlas, even though some of the methods we consider extract the regions from the data rather than relying on a reference atlas (see Figure 1). We investigate here pipelines that discriminate individuals based on the connection strength of this connectome [START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF], with a classifier on the edge weights [START_REF] Craddock | Disease state prediction from resting state functional connectivity[END_REF].

Specifically, we consider connectome-extraction pipelines composed of four steps: 1) region estimation, 2) time series extraction, 3) matrix estimation, and 4) classification -see Figure 1. We investigate different options for the various steps: brain parcellation method, connectivity matrix estimation, and final classifier.

In order to measure the impact of uncontrolled variability on prediction performance, we designed a diagnosis task where the model classifies participants coming from an acquisition site not seen during training.

Datasets: Preprocessing and Quality Assurance (QA)

We ran our analyses on the ABIDE repository [START_REF] Martino | The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism[END_REF], a sample aggregated across 17 independent sites. Since there was no prior coordination between sites, the scan and diagnostic/assessment protocols vary across sites1 . To easily replicate and extend our work, we rely on a publicly available preprocessed version of this dataset provided by the Preprocessed Connectome Project2 initiative. We specifically used the data processed with the Configurable Pipeline for the Analysis of Connectomes [START_REF] Craddock | Towards automated analysis of connectomes: The configurable pipeline for the analysis of connectomes (c-pac)[END_REF] (C-PAC). Data were selected based on the results of quality visual inspection by three human experts who inspected for largely incomplete brain coverage, high movement peaks, ghosting and other scanner artifacts. This yielded 871 subjects out of the initial 1112.

Preprocessing included slice-timing correction, image realignment to correct for motion, and intensity normalization. Nuisance regression was performed to remove signal fluctuations induced by head motion, respiration, cardiac pulsation, and scanner drift [START_REF] Lund | Motion or activity: their role in intra-and inter-subject variation in fMRI[END_REF][START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF]. Head motion was modeled using 24 regressors derived from the parameters estimated during motion realignment [START_REF] Friston | Statistical parametric maps in functional imaging: A general linear approach[END_REF], scanner drift was modeled using a quadratic and linear term, and physiological noise was modeled using the 5 principal components with highest variance from a decomposition of white matter and CSF voxel time series (CompCor) [START_REF] Behzadi | A component based noise correction method (CompCor) for bold and perfusion based fMRI[END_REF]. After regression of nuisance signals, fMRI was coregistered on the anatomical image with FSL BBreg, and the results where normalized to MNI space with the non-linear registration from ANTS [START_REF] Avants | Advanced normalization tools (ants)[END_REF]. Following time series extraction, data were detrended and standardized (dividing by the standard deviation in each voxel).

Step 1: Region definition

To reduce the dimensionality of the estimated connectomes, and to improve the interpretability of results, connectome nodes were defined from regions of interest (ROIs) rather than single voxels. ROIs can be either selected from a priori atlases or directly estimated from the correlation structure of the data being analyzed. Several choices for both approaches were used and compared to evaluate the impact of ROI definition on prediction accuracy. Note that the ROI were all defined on the training set to avoid possible overfitting.

We considered the following predefined atlases: Harvard Oxford [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF], a structural atlas computed from 36 individuals' T1 images, Yeo [START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF], a functional atlas composed of 17 networks derived from clustering [START_REF] Lashkari | Discovering structure in the space of fMRI selectivity profiles[END_REF] of functional connectivity on a thousand individuals, and Craddock [START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF], a multiscale functional atlas computed using constrained spectral clustering, to study the impact of the number of regions.

We derived data-driven atlases based on four strategies. We explored two clustering methods: K-Means, a technique to cluster fMRI time series [START_REF] Goutte | On clustering fMRI time series[END_REF][START_REF] Thirion | Which fMRI clustering gives good brain parcellations?[END_REF], which is a top-down approach minimizing cluster-level variance; and Ward's clustering, which also minimizes a variance criterion using hierarchical agglomeration. Ward's algorithm admits spatial constraints at no cost and has been extensively used to learn brain parcellations [START_REF] Thirion | Which fMRI clustering gives good brain parcellations?[END_REF]. We also explored two decomposition methods: Independent component analysis (ICA) and multi-subject dictionary learning (MSDL). ICA is a widely-used method for extracting brain maps from R-fMRI [START_REF] Calhoun | A method for making group inferences from fMRI data using independent component analysis[END_REF][START_REF] Beckmann | Probabilistic independent component analysis for functional magnetic resonance imaging[END_REF] that maximizes the statistical independence between spatial maps. We specifically used the ICA implementation of [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF]. MSDL extracts a group atlas based on dictionary learning and multi-subject modeling, and employs spatial penalization to promote contiguous regions [START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF]. Unlike MSDL and Ward clustering, ICA and K-Means do not take into account the spatial structure of the features and are not able to enforce a particular spatial constraint on their components. We indirectly provide such structure by applying a Gaussian smoothing on the data [START_REF] Thirion | Which fMRI clustering gives good brain parcellations?[END_REF] with a full-width half maximum varying from 4mm to 12mm. Differences in the number of ROIs in the various atlases present a potential confound for our comparisons. As such, connectomes were constructed using only the 84 largest ROIs from each atlas, following [49] 3 .

Step 2: Time-series extraction

The common procedure to extract one representative time series per ROI is to take the average of the voxel signals in each region (for non overlapping ROIs) or the weighted average of non-overlapping regions (for fuzzy maps). For example, ICA and MSDL produce overlapping fuzzy maps. Each voxel's intensity reflects the level of confidence that this voxel belongs to a specific region in the atlas. In this case, we found that extracting time courses using a spatial regression that models each volume of R-fMRI data as a mixture of spatial patterns gives better results (details in Appendix A and Figure A3). For nonoverlapping ROIs, this procedure is equivalent to weighted averaging on the ROI maps.

After this extraction, we remove at the region level the same nuisance regressors as at the voxel level in the preprocessing. Signals summarizing high-variance voxels (CompCor [START_REF] Behzadi | A component based noise correction method (CompCor) for bold and perfusion based fMRI[END_REF]) are regressed out: the 5 first principal components of the 2% voxels with highest variance. In addition, we extract the signals of ROIs corresponding to physiological or acquisition noise 4 and regress them out to reject non-neural information [START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF]. Finally, we also regress out primary and secondary head motion artifacts using 24regressors derived from the parameters estimated during motion realignment [START_REF] Friston | Statistical parametric maps in functional imaging: A general linear approach[END_REF].

Step 3: Connectivity matrix estimation

Since scans in the ABIDE dataset do not include enough R-fMRI time points to reliably estimate the true 3 For Harvard-Oxford and Yeo reference atlases, selecting 84 ROIs yields slightly better results than using the full version as shown in Figure A2. This is probably due to the fact that these additional regions are smaller and less important regions, hence they induce spurious correlations in the connectivity matrices. 4 We identify these ROIs using an approach similar to FSL FIX (FMRIB's ICA-based Xnoiseifier). FIX extracts spatial and temporal descriptors from matrix decomposition results and uses a classifier trained on manually labeled components to determine if they correspond to noise.

covariance matrix for a given individual, we used the Ledoit-Wolf shrinkage estimator, a parameter-free regularized covariance estimator [START_REF] Ledoit | A well-conditioned estimator for largedimensional covariance matrices[END_REF], to estimate relationships between ROI time series (details are given in Appendix B). For pipeline step 3, we studied three different connectivity measures to capture interactions between brain regions: i) the correlation itself, ii) partial correlations from the inverse covariance matrix [START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF][START_REF] Smith | Network modelling methods for fMRI[END_REF], and iii) the tangent embedding parameterization of the covariance matrix [START_REF] Varoquaux | A group model for stable multi-subject ICA on fMRI datasets[END_REF][START_REF] Ng | Transport on riemannian manifold for functional connectivity-based classification[END_REF]. We obtained one connectivity weight per pair of regions for each subject. At the group level, we ran a nuisance regression across the connectivity coefficients to remove the effects of site, age and gender.

Step 4: Supervised learning

As ABIDE represents a post-hoc aggregation of data from several different sites, the assessments used to measure autism severity vary somewhat across sites. As a consequence, autism severity scores are not directly quantitatively comparable between sites. However, the diagnosis between ASD and Typical Control (TC) is more reliable and has already been used in the context of classification [START_REF] Nielsen | Multisite functional connectivity MRI classification of autism: ABIDE results[END_REF][START_REF] Plitt | Functional connectivity classification of autism identifies highly predictive brain features but falls short of biomarker standards[END_REF][START_REF] Chen | Diagnostic classification of intrinsic functional connectivity highlights somatosensory, default mode, and visual regions in autism[END_REF].

Discriminating between ASD and TC individuals is a supervised-learning task. We use the connectivity measures between all pairs of regions, extracted in the previous step, as features to train a classifier across individuals to discriminate ASD patients from TC.

We explore different machine-learning methods 5 . We first rely on the most commonly used approach, the 2penalized support vector classification SVC. Given that we expect only a few connections to be important for diagnosis, we also include the 1 -penalized sparse SVC. Finally, we include ridge regression, which also uses an 2 penality but is faster than the SVC. For all models we use the implementation of the scikit-learn library [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] 6 . We chose to work only with linear models for their interpretability. In particular, in subsection 3.7 we inspect classifier weights to identify which connections are the most important for prediction.

Exploring dataset heterogeneity

Previous studies of prediction using the ABIDE sample [START_REF] Plitt | Functional connectivity classification of autism identifies highly predictive brain features but falls short of biomarker standards[END_REF][START_REF] Chen | Diagnostic classification of intrinsic functional connectivity highlights somatosensory, default mode, and visual regions in autism[END_REF] have included less than 25% of the available data, most likely to limit heterogeneity -not only due to differences in imaging, but also to factors such as age, sex, and handedness. We explored the effect of such choice on the results, by examining data subsets of different heterogeneity (see Table 1). These included: subsample #1 - A2 for more details) 7 .

2.7.

Experiments: prediction on heterogeneous data ABIDE datasets come from 17 international sites with no prior coordination, which is a typical clinical application setting. In this setting, we want to: i) measure the performance of different prediction pipelines, ii) use these measures to find the best options for each pipeline step (Figure 1) and, iii) extract ASD biomarkers from the best pipeline.

Cross-validation. In order to measure the prediction accuracy of a pipeline, we use 10-fold cross-validation by training it exclusively on training data (including atlas estimation) and predicting on the testing set. Any hyperparameter of the various methods used inside a pipeline is set internally in a nested cross-validation.

We used two validation approaches. Inter-site prediction addresses the challenges associated with aggregate datasets by using whole sites as testing sets. This scenario also closely simulates the real world clinical challenge of not being able to sample data from every possible site. Note that this cross-validation can only be applied on a dataset with at least 10 acquisition sites, in order to leave out a different site in each fold. We do not apply it on subsample #5 that has been restricted to 3 sites to reduce sitespecific variability. Intra-site prediction builds training and testing sets are as homogeneous as possible, reducing site-related variability. It is based on stratified shuffle split cross-validation, which splits participants into training and test sets while preserving the ratio of samples for each site and condition. We used 80% of the dataset for training and the remaining for testing.

Learning curves. A learning curve measures the impact of the number of training samples on the performance of a predictive model. For a given prediction task, the learning curve is constructed by varying the number of samples in the training set with a constant test set. Typically, an increasing curve indicates that the model can still gain additional useful information from more data. It is expected that the performance of a model will plateau at some point.

Summarizing prediction scores for a pipeline choice. To quantify the impact of the different options on the prediction accuracy, we select representative prediction scores for each model. Indeed, for ICA, K-Means, and MSDL, we explore several values for each parameters and thus obtain a range of prediction values. For these methods, we use the 10% best-performing pipelines for post-hoc analysis8 . We do not rely only on the top-performing pipelines, as they may result from overfitting, nor the complete set of pipelines as some pipelines yield bad scores because of bad parametrization and may artificially lower the performance of these methods.

Impact of pipeline choices on prediction accuracy. In order to determine the best pipeline option at each step, we want to measure the impact of these choices on the prediction scores. In a full-factorial ANOVA, we then fit a linear model to explain prediction accuracy from the choice of pipeline steps, modeled as categorical variables. We report the corresponding effect size and its 95% confidence interval.

Pairwise comparison of pipeline choice. Two options can also be compared by comparing the scores between pairs of pipelines that differ only by that option using a Wilcoxon test that does not rely on Gaussian assumptions.

Extracting biomarkers. For a given pipeline, the final classifier yields a biomarker to diagnose between HC and ASD. To characterize this biomarker, we measure p-values on the classifier weights. We obtain the null distribution of the weights by permuting the prediction labels.

Computational aspects

For data-driven brain atlas models (e. g. ICA), regionextraction methods entail feature learning on up to 300 GB of data, too big to fit in memory. Systematic study requires many cross-validations and subsamples and thus careful computational optimizations. Hence, for K-Means and ICA, we reduced the length of the time series by applying PCA dimensionality reduction with efficient randomized SVD [START_REF] Martinsson | A randomized algorithm for the decomposition of matrices[END_REF]. MSDL uses algorithmic optimizations [START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF] to process a large number of subjects in a reasonable time 9 (2 hours for 871 subjects). However, it is still the most costly method because of the number of runs required to find the optimal value for its 3 parameters.

For reproducibility, we rely on the publicly available implementations of the scikit-learn package [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] for all the clustering algorithms, connectivity matrix estimators and predictors.

Results

Here we present the most notable trends emerging from our analyses. More details are provided in supplementary materials. First, we compared inter-and intra-site prediction of ASD diagnostic labels while varying the number of subjects in the training set. Second, in a post-hoc analysis, we identified the pipeline choices most relevant for prediction and proposed a good choice of pipeline steps for prediction. Finally, by analyzing the weights of the classifiers, we highlighted functional connections that best discriminate between ASD and TC.

Overall prediction results

For inter-site prediction, the highest accuracy obtained on the whole dataset is 66.8% (see Table 2) which exceeds previously published ABIDE findings [START_REF] Nielsen | Multisite functional connectivity MRI classification of autism: ABIDE results[END_REF][START_REF] Plitt | Functional connectivity classification of autism identifies highly predictive brain features but falls short of biomarker standards[END_REF][START_REF] Chen | Diagnostic classification of intrinsic functional connectivity highlights somatosensory, default mode, and visual regions in autism[END_REF]] and chance at 53.7% 10 . Importantly, performance increases steadily with training set size, for all subsamples (Figures 2 andA6). This increase implies that optimal classification performance is not reached even for the largest dataset tested.

The main difference between intra-site and inter-site settings is that the variability of inter-site prediction is higher (Table 2). However, this difference disappears with a large enough training set (Figure 2). This is encouraging for clinical applications.

As shown in Figure 3, the best predicting pipeline is MSDL, tangent space embedding and l 2 -regularized classifiers (SVC-l 2 and ridge classifier). All effects are observed 9 Despite the optimizations, the whole study presented here is very computationally costly. Indeed, the nested cross-validation, the various pipeline options, and the various subsets, would lead to a computation time of 10 years on a computer with 8 cores and 32GB RAM. We used a computer grid to reduce it to two months. 10 Chance is calculated by taking the highest score obtained by a dummy classifier that predicts random labels following the distribution of the two classes. with p < 0.01. Detailed pairwise comparisons (Figure 4) confirm these trends by showing the effect of each option compared to the best ones. Results of intra-site prediction have smaller standard error, confirming higher variability of inter-site prediction results (Figure A5).

Subsample

Effect of the choice of atlas

The choice of atlas appears to have the greatest impact on prediction accuracy (Figure 3). Estimation of an atlas from the data with MSDL, or using reference atlases leads to maximal performance. Of note, for smaller samples, data-driven atlas-extraction methods other than MSDL perform poorly (Figure A8). In this regime, noise

Cross validation scheme ABIDE subsample

Step 1: Atlas extraction

Step 3: Covariance matrix

Step 4: Predictor We also prepended two steps corresponding to the cross-validation schemes and ABIDE subsamples. Each bar represents the impact of the corresponding option on the prediction accuracy, relative to the mean prediction. This effect is measured via a full-factorial analysis of variance (ANOVA), analyzing the contribution of each step in a linear model. Each step of the pipeline is considered as a categorical variable. Error bars give the 95% confidence interval. Multi Subject Dictionary Learning (MSDL) atlas extraction method gives significantly better results while reference atlases are slightly better than the mean. Among all matrix types, tangent embedding is the best on all ABIDE subsets. Finally, l 2 -regularized classifiers perform better than l 1 -regularized ones. Results for each subsamples are presented in the appendices -Figure A8. can be problematic and limit the generalizability of atlases derived from data.

MSDL performs well regardless of sample size, and the performance gain is significant compared to all other strategies apart from the Harvard-Oxford atlas (Figure 4). This is likely attributable to MSDL's strong spatiallystructured regularization that increases its robustness to noise.

Effect of the covariance matrix estimator

While correlation and partial correlation approaches are currently the two dominant approaches to connectome estimation in the current R-fMRI literature, our results highlight the value of tangent-space projection that outperforms the other approaches in all settings (Figure 3 and Figure 4), though the difference is significant only compared with correlation matrices (Figure 4). This gain in prediction accuracy is not without a cost, as the matrices derived cannot be read as correlation matrices. Yet, the differences in connectivity that they capture can be directly interpreted [START_REF] Varoquaux | Detection of brain functional-connectivity difference in post-stroke patients using group-level covariance modeling[END_REF].

Among the correlation-based approaches, full correlations outperformed partial correlations, most notably for intra-site prediction. These limitations of partial correlations may reflect estimation challenges 11 with the relatively short time series included in ABIDE datasets (typically 5-6 minutes per participant).

Effect of the final classifier

The results show that 2 -regularized classifiers perform best in all settings (Figure 3 and Figure 4). This may be due to the fact that a global hypoconnectivity, as previously observed in ASD patients, is not well captured by

Effect of the dataset size and heterogeneity

While not a property of the data-processing pipeline itself, a large number of training subjects is the most important factor of success for prediction. Indeed, we observe that prediction results improve with the number of subjects, even for large number of subjects (Figure 3). Despite that increase, we note that prediction on the largest subsample (subsample #1) gives lower scores than prediction on subsample #2, which may be due to the variance introduced by small sites. We also note that the importance of the choice between each option of the pipeline decreases with the number of subjects (Figure A8).

Number of regions

Effect of the number of regions

In order to compare models of similar complexity, the previous experiments were restricted to 84 regions as in [START_REF] Abraham | Extracting brain regions from rest fMRI with Total-Variation constrained dictionary learning[END_REF]. However, this choice is based on previous studies and may not be the optimal choice for our problem. To validate this aspect of the study, we also varied the number of regions and analyzed its impact on the classification results in Figure 5, similarly to Figure 3. In order to avoid computational cost 12 , we explored only two atlases: the data-driven Ward's hierarchical clustering and spectral clustering atlases computed in [START_REF] Craddock | A whole brain fMRI atlas generated via spatially constrained spectral clustering[END_REF].

Across all settings, we distinguish 3 different regimes. First, below 20 ROIs, performance is bad. Indeed very large regions average out signals of different functional areas. Above 140 ROIs, the results seem to be unstable, though this trend is less clear in intra-site prediction. Finally, the best results are between 40 and 100 ROIs. Our choice of 84 ROIs is within the range of best performing models. 12 As MSDL is much more computationally expensive than Ward, we have not explored the effect of modifying of the number of regions extracted with MSDL. However, we expect that the optimal number of regions is not finely dependent on the region-extraction method.

Inspecting discriminative connections

A connection between two regions is considered discriminative if its value helps diagnosing ASD from TC. To understand what is driving the best prediction pipeline, we extract its most discriminative connections as described in section 2.7. Given that the 10-fold cross-validation yields 10 different atlases, we start by building a consensus atlas by taking all the regions with a DICE similarity score above 0.9. We obtain an atlas composed of 37 regions (see Figure 6). These regions and the discriminative connections form the neurophenotypes. Note that discriminative connections should be interpreted with caution. Indeed, our study covers a wide range of age and diagnostic criteria. Additionally, predictive models cannot be interpreted as edge-level tests.

Default Mode Network (Figure 6.a). We observe a lower connectivity between left and right temporo-parietal junctions. Both hypoconnectivity between regions on adolescent and adult patients [START_REF] Cherkassky | Functional connectivity in a baseline resting-state network in autism[END_REF][START_REF] Kennedy | Failing to deactivate: resting functional abnormalities in autism[END_REF] and hyperconnectivity in children [START_REF] Supekar | Brain hyperconnectivity in children with autism and its links to social deficits[END_REF] have been reported. Our findings are concordant with [START_REF] Monk | Abnormalities of intrinsic functional connectivity in autism spectrum disorders[END_REF] that observed lower fronto-parietal connectivity and stronger connectivity between temporal ROIs.

Pareto-insular network (Figure 6.b). We observe interhemispheric hypoconnectivity between left and right anterior insulae, and left and right inferior parietal lobes. Those regions are part of a frontoparietal control network related to learning and memory [START_REF] Iidaka | Frontoparietal network involved in successful retrieval from episodic memory. spatial and temporal analyses using fMRI and ERP[END_REF]. Such hypoconnectivity has already been widely evoked in ASD [START_REF] Anderson | Decreased interhemispheric functional connectivity in autism[END_REF][START_REF] Martino | Aberrant striatal functional connectivity in children with autism[END_REF] and previously found in the ABIDE dataset [START_REF] Plitt | Functional connectivity classification of autism identifies highly predictive brain features but falls short of biomarker standards[END_REF][START_REF] Ray | Structural and functional connectivity of the human brain in autism spectrum disorders and attentiondeficit/hyperactivity disorder: A rich club-organization study[END_REF]. These ROIs also match anatomical regions of increased cortical density in ASD [START_REF] Haar | Anatomical abnormalities in autism?[END_REF].

Semantic ROIs (Figure 6.c). We observe a more complicated pattern in ROIs related to language. Connectivity is stronger between the right supramarginal gyrus and the left Middle Temporal Gyrus (MTG) while a lower connectivity is observed between the right MTG and the left temporo-parietal junction. Differences of laterality in language networks implicating the middle temporal gyri have already been observed [START_REF] Kleinhans | Atypical functional lateralization of language in autism spectrum disorders[END_REF]. Although this symptom is often considered as decorrelated from ASD, we find that these regions are relevant for diagnosis.

Discussion

We studied pipelines that extract neurophenotypes from aggregate R-fMRI datasets through the following steps: 1) region definition from R-fMRI, 2) extraction of regions activity time series 3) estimation of functional interactions between regions, and 4) construction of a discriminant model for brain-based diagnostic classification. The proposed pipelines can be built with the Nilearn neuroimaging-analysis software and the atlas computed with MSDL is available for download 13 . We have demonstrated that they can successfully predict diagnostic status across new acquisition sites in a real-world situation, a large Autism database. This validation with a leave-onesite-out cross-validation strategy reduces the risk of overfitting data from a single site, as opposed to the commonly used leave-one-subject-out validation approach [START_REF] Nielsen | Multisite functional connectivity MRI classification of autism: ABIDE results[END_REF]. It enabled us to measure the impact of methodological choices on prediction. In the ABIDE dataset, this approach classified the label of autism versus control with an accuracy of 67% -a rate superior to prior work using the larger ABIDE sample. The heterogeneity of weakly-controlled datasets formed by aggregation of multiple sites poses great challenges to develop brain-based classifiers for psychiatric illnesses. Among those most commonly cited are: i) the many sources of uncontrolled variation that can arise across studies and sites (e. g. scanner type, pulse sequence, recruitment strategies, sample composition), ii) the potential for developing classifiers that are reflective of only the largest sites included, and iii) the danger of extracting biomarkers that will be useful only on sites included in the training sample. Counter to these concerns, our results demonstrated the feasibility of using weakly-controlled heterogeneous datasets to identify imaging biomarkers that are robust to differences among sites. Indeed, our prediction accuracy of 67% is the highest reported to-date on the whole ABIDE dataset (compared to 60% in [START_REF] Nielsen | Multisite functional connectivity MRI classification of autism: ABIDE results[END_REF]), though 13 http://team.inria.fr/parietal/files/2015/07/MSDL_ABIDE. zip better prediction has been reported on curated more homogeneous subsets [START_REF] Plitt | Functional connectivity classification of autism identifies highly predictive brain features but falls short of biomarker standards[END_REF][START_REF] Chen | Diagnostic classification of intrinsic functional connectivity highlights somatosensory, default mode, and visual regions in autism[END_REF]. Importantly, we found that increasing sample size helped tackling heterogeneity as intersite prediction accuracy approached that of intra-site prediction. This undoubtedly increases the perceived value that can be derived from aggregate datasets, such as the ABIDE and the ADHD-200 [START_REF] Consortium | The adhd-200 consortium: a model to advance the translational potential of neuroimaging in clinical neuroscience[END_REF].

Default Mode Network L R L R a Pareto-insular network L R L R b Semantic ROIs L R L R c L R y=-54 x=0 L R z=3
Our results suggest directions to improve prediction accuracy. First, more data is needed to optimize the classifier. Indeed, our experiments show that with a training set of 690 participants aggregated across multiple sites (80% of 871 available subjects), the pipelines have not reached their optimal performance. Second, a broader diversity of data will be needed to more definitively assess prediction for sites not included in the training set. Both of these needs motivate more data sharing.

Methodologically, our systematic comparison of choices in the pipeline steps can give indications to establish standard processing pipelines. Given that processing pipelines differ greatly between studies [START_REF] Carp | The secret lives of experiments: methods reporting in the fMRI literature[END_REF], comparing data-analysis methods is challenging. On the ABIDE dataset, we found that a good choice of processing steps can strongly improve prediction accuracy. In particular, the choice of regions is most important. We found that extracting these regions with MSDL gives best results. However, the benefit of this method is not always clear cut -on small datasets, references atlases are also good performers. For the rest of the pipeline, tangent embedding and 2 -regularized classifier are overall the best choice to maximize prediction accuracy.

While the primary goals of this study were methodo-logical, it is worth noting that biomarkers identified were largely consistent with the current understanding of the neural basis of ASD. The most informative features to predict ASD concentrated in the intrinsic functional connectivity of three main functional systems: the defaultmode, parieto-insular, and language networks, previously involved in ASD [START_REF] Vissers | Brain connectivity and high functioning autism: a promising path of research that needs refined models, methodological convergence, and stronger behavioral links[END_REF][START_REF] Rane | Connectivity in autism: A review of mri connectivity studies[END_REF][START_REF] Hernandez | Neural signatures of autism spectrum disorders: insights into brain network dynamics[END_REF]. Specifically, decreased homotopic connectivity between temporo-parietal junction, insula and inferior parietal lobes appeared to characterize ASD. Decreases in homotopic functional connectivity have been previously reported in R-fMRI studies with moderately [START_REF] Anderson | Decreased interhemispheric functional connectivity in autism[END_REF][START_REF] Dinstein | Disrupted neural synchronization in toddlers with autism[END_REF] or large samples such as ABIDE [START_REF] Martino | The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism[END_REF].

Here our findings point toward a set of inter-hemispheric connections involving heteromodal association cortex subserving social cognition (i. e. temporo-parietal junction [START_REF] Saxe | People thinking about thinking people: the role of the temporo-parietal junction in theory of mind[END_REF]), and cognitive control (anterior insula and right inferior parietal cortex [START_REF] Menon | Saliency, switching, attention and control: a network model of insula function[END_REF]), commonly affected in ASD [START_REF] Anderson | Decreased interhemispheric functional connectivity in autism[END_REF][START_REF] Martino | Aberrant striatal functional connectivity in children with autism[END_REF][START_REF] Plitt | Functional connectivity classification of autism identifies highly predictive brain features but falls short of biomarker standards[END_REF][START_REF] Ray | Structural and functional connectivity of the human brain in autism spectrum disorders and attentiondeficit/hyperactivity disorder: A rich club-organization study[END_REF]. Limitations of our biomarkers may arise from the heterogeneity of ASD as a neuropsychiatric disorder [START_REF] Geschwind | Autism spectrum disorders: developmental disconnection syndromes[END_REF].

Beyond forming a neural correlate of the disorder, predictive biomarkers can help define diagnosis subgroupsthe next logical endeavor in classification studies [START_REF] Loth | Identification and validation of biomarkers for autism spectrum disorders[END_REF].

As a psychiatric study, the present work has a number of limitations. First, although composed of samples from geographically distinct sites, the representativeness of the ABIDE sample has not been established. As such, the features identified may be biased and necessitate further replication; though their consistency with prior findings alleviates this concern to some degree. Further work calls for characterizing this prediction pipeline on more participants (e. g. ABIDE II) and different pathologies (e. g. ADHD-200 dataset). Second, while the present work primarily relied on prediction accuracy to assess the classifiers derived, a variety of alternative measures exist (e. g. sensitivity, specificity, J-statistic). This decision was largely motivated to facilitate comparison of our findings with those from prior studies using ABIDE, which most frequently used prediction accuracy. Depending on the specific application for a classifier, the prediction metric to optimize may vary. For example screening tools would need to favor sensitivity over accuracy or specificity; in contrast with medical diagnostic tests that require specificity over sensitivity or accuracy [START_REF] Castellanos | Clinical applications of the functional connectome[END_REF]. The versatility of the pipeline studied allows us to maximize either of these scores depending on the experimenter's needs. Finally, while the present work focused on the ABIDE dataset, it is likely that its findings regarding optimal pipeline decisions will carry over to other aggregate samples, as well as more homogeneous samples. With that said, future applications will be required to verify this point.

In sum, we experimented with a large number of different pipelines and parametrizations on the task of predicting ASD diagnosis on the ABIDE dataset. From an extensive study of the results, we have found that R-fMRI from a large amount of participants sampled from many differ-ent sites could lead to functional-connectivity biomarkers of ASD that are robust, including to inter-site variability. Their predictive power markedly increases with the number of participants included. This increase holds promise for optimal biomarkers forming good probes of disorders.

With w ∈ R n , x ∈ R n×m , and y ∈ R m . In our case, n, the number of features, is the number of connections between pairs of brain regions and m, the number of samples, is the number of subjects. We explore the sparse regularization 1 and 2 -penalized SVC (resp. p = 1 and p = 2). We use the implementation of scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] that relies on the implementation provided by LibLinear [START_REF] Fan | LIBLINEAR: A library for large linear classification[END_REF].

Ridge classifier

Ridge regression can be used in the context of classification by attributing binary values to each of the classes (typically -1 and 1). Classifying an sample boils down to looking at the sign of the prediction of the regressor.

The Ridge regression minimizes the square distance to the hyperplane, like an ordinary least squares, but it adds an 2 penalization on the weights in order to reward coefficients close to 0 and thus obtain the smallest possible model. This prevents colinear feature to skew the coefficients of the model. The minimized function is:

min w y -w T x 2 + w 2
We use the implementation of scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF].

Appendix D. Effect of movements

Participant's movements during fMRI acquisition is known to induce spurious correlations in the connectome estimation. In our particular case, we expect healthy individuals and ASD ones to exhibit different movement patterns because of the behavioral symptoms of the latter.

In our study, we use standard methods to remove the effect of movement as much as possible. We also performed an analysis to measure the amount of information contained in the movement patterns and ensure that our predictions are not based on potential movement residuals.

Movement regression

In the experiments, during time series extraction (step 2), we regress out all possible confounding factors, which includes movement estimation. In order to quantify the information contained in the movements, we try predicting ASD diagnosis with and without movement regression. Figure A1 shows that no significant difference is observed with and without movement regression, which means that movement effect does not affect strongly our connectome estimation.

Movement-based prediction

Going further, we also perform the diagnosis directly based on movement estimations. Because movement parameters are time-series, it is not possible to use them directly as features in a prediction task. For this purpose, we extract temporal descriptors from these series. These descriptors are similar to the ones extracted in FSL FIX [START_REF] Salimi-Khorshidi | Automatic denoising of functional MRI data: combining independent component analysis and hierarchical fusion of classifiers[END_REF], a tool to identify noisy ICA components: coefficients of autoregressive models of several orders, kurtosis, skewness, entropy, difference between the mean and the median, and Fourier transform coefficients.

As for the prediction pipelines, we tried a wide range of strategies to predict from movement parameters. We extracted the descriptors on the movements, the gradient of the movements, the gradient squared, and the framewise displacement. We tried several predictors (with and without feature extraction), logistic regression, SVC, Ridge classifier, Gaussian Naive Bayes and random forests.

We present the results obtained on the best prediction procedure. For each participant, 56 descriptors are extracted from the movement estimations and their gradient, and are used as features in an SVC for diagnosis. We observe that the scores obtained using movements are at chance level. Note that we see a trend already observed on prediction based on connectomes: The variance of prediction scores is much larger for inter-site prediction. A5: Average sensitivity scores (along with their standard deviations) for top performing pipelines. This table summarizes the scores for pipelines with best prediction accuracy for each atlas and subset using intra-site or inter-site prediction. Best results are shown in bold.

Intra-site Inter-site

Intra-site Inter-site Intra-site Inter-site

Figure 1 :

 1 Figure1: Functional MRI analysis pipeline. Cross-validation schemes used to validate the pipeline are presented above. Intra-site cross-validation consists of randomly splitting the participants into training and testing sets while preserving the ratio of samples for each site and condition. Inter-site cross-validation consists of leaving out participants from an entire site as testing set. In the first step of the pipeline, regions of interest are estimated from the training set. The second step consists of extracting signals of interest from all the participants, which are turned into connectivity features via covariance estimation at the third step. These features are used in the fourth step to perform a supervised learning task and yield an accuracy score. An example of pipeline is highlighted in red. This pipeline is the one that gives best results for inter-site prediction. Each model is decribed in the section relative to material and methods.

Figure 2 :

 2 Figure2: Learning curve. Classification results obtained by varying the ratio of the training set using to train the classifier while keeping a fixed testing set. The colored bands represent the standard error of the prediction. A score increasing with the number of subjects, i. e. a positive slope indicates that the addition of new subjects improves performance. This curve is an average of the results obtained on several subsamples. Detailed results per subsamples are presented in the appendices -FigureA6

Figure 3 :

 3 Figure 3: Impact of pipeline steps on prediction. Each block of bars represents a step of the pipeline (namely step 1, 3 and 4).We also prepended two steps corresponding to the cross-validation schemes and ABIDE subsamples. Each bar represents the impact of the corresponding option on the prediction accuracy, relative to the mean prediction. This effect is measured via a full-factorial analysis of variance (ANOVA), analyzing the contribution of each step in a linear model. Each step of the pipeline is considered as a categorical variable. Error bars give the 95% confidence interval. Multi Subject Dictionary Learning (MSDL) atlas extraction method gives significantly better results while reference atlases are slightly better than the mean. Among all matrix types, tangent embedding is the best on all ABIDE subsets. Finally, l 2 -regularized classifiers perform better than l 1 -regularized ones. Results for each subsamples are presented in the appendices -FigureA8.

Figure 4 :

 4 Figure 4: Comparison of pipeline options. Each plot compares the classification accuracy scores if one parameter is changed and the others kept fixed. We measured statistical signficance using a Wilcoxon signed rank test and corrected for multiple comparisons. Detailed one-to-one comparison of the scores are shown in the appendices -Figure A5

Figure 5 :

 5 Figure 5: Impact of the number of regions on prediction: each bar indicates the impact of the number of regions on the prediction accuracy relative to the mean of the prediction. These values are coefficients in a linear model explaining the best classification scores as function of the number of regions. Error bars give the 95% confidence interval, computed by a full factorial ANOVA. Atlases containing more than 40 ROIs give better results in all settings. Results per subsamples are presented in the appendices -Figure A9.

Figure 6 :

 6 Figure 6: Significant non-zero connections in the predictive biomarkers distinguishing controls from ASD patients. Red connections are stronger in controls and blue connections are stronger in ASD patients. Subfigures a, b, and c reported intranetwork difference. On the left is the consensus atlas obtained by selecting regions consistently extracted on 10 subsets of ABIDE. Colors are randomly assigned. Results obtained on other networks are shown in the appendices -Figure A7.

Figure A5 :

 A5 Figure A5: Comparison of pipeline options. Each plot compares the results of the best option against the other for each step. Points above the line means that the option show on y-axis is better than option on the x-axis. The 2 cross validation schemes are represented by circles and stars.

Figure A6 :

 A6 Figure A6: Learning curves. Classification results obtained by varying the number of subjects in the training set while keeping a fixed testing set. The colored band represent the standard error of the prediction accuracy. A score increasing with the number of subjects, i. e. an increasing curve, indicates that the addition of new subjects brings information.

  Figure A7: Connections significantly non-zero in the predictive biomarkers distinguishing controls from ASD patients. Red connections are stronger in controls and blue connections are stronger in ASD patients. Subfigures a, b, c, d and e reported intra-network difference. below is the consensus atlas extracted by selecting regions consistently extracted on 10 subsets of ABIDE. Colors are random.L R

Table 1 :

 1 Subsets of ABIDE used in evaluation. We explore several subsets of ABIDE with different homogeneity. Card. stands for the cardinality of the dataset, i. e. the number of subjects. More details about these subsets are presented in the appendices -TableA2.all participants (871 individuals, 727 males and 144 females, 6 to 64 years old) is the full set of individuals that passed QA, subsample #2 -largest sites (736 individ-

	Subsample		Card. Sites	Selection criteria
	1 All subjects	871	16	All subjects after QC
	2 Biggest sites	736	11	Remove 6 smallest sites
	3 Right	handed	639	14	All subjects without left
	males				handed and women
	4 Right	handed	420	14	Right handed males be-
	males, 9-18 yo			tween 9 and 18 years old
	5 Right	handed	226	3	Young right handed
	males, 9-18 yo, 3			males from 3 biggest
	sites				sites

uals, 613 males and 123 females) 6 sites with less than 30 participants are removed from subsample #1, subsample #3 -right handed males (639 individuals) consists of the right-handed males from subsample #1, subsample #4 -right handed males, 9-18 years (420 individuals) is the restriction of subsample #3 to participants between 9 and 18 years old; subsample #5 -right handed males, 9-18 yo, 3 sites (226 individuals) are the individuals from subsample #4 belonging to the three largest sites (see Table

Table 2 :

 2 Average accuracy, specificity and sensitivity scores (and standard deviation) for top performing pipelines. Accuracy is the fraction of well-classified individuals (chance level is at 53.7%). Specificity is the fraction of well classified healthy individuals among all healthy individuals (percentage of well-classified negatives). Sensitivity is the ratio of ASD individuals among all subjects classified as ASD (percentage of positives that are true positives). Results per atlas are available in the appendices -TableA5.

		#5	#4	#3	#2	#1
	Inter-site		69.7%	65.1%	68.7%	66.8%
	Accuracy		±8.9% ±5.8% ±9.3% ±5.4%
	Inter-site		74.0%	69.1%	73.9%	72.3%
	Specificity		±12.9%	±7.8%	±11.7% ±12.1%
	Inter-site		65.4%	61.3%	62.8%	61.0%
	Sensitivity		±13.7%	±7.7%	±13.6% ±17.0%
	Intra-site	66.6%	65.8%	65.7%	67.9%	66.9%
	Accuracy	±5.4% ±5.9% ±4.9% ±1.9% ±2.7%
	Intra-site	78.4%	72.3%	77.5%	76.6%	78.3%
	Specificity	±9.7%	±9.1%	±8.1%	±5.0%	±4.1%
	Intra-site	53.2%	58.1%	52.2%	58.1%	53.2%
	Sensitivity	±10.9%	±9.2%	±9.1%	±5.0%	±5.8%

Table A1 :

 A1 Prediction accuracy scores based on movements estimations. Results are similar across folds but variance of intersite prediction is larger. Comparison of prediction results obtained using 84 regions extracted from Harvard Oxford and Yeo atlases the set of all regions extracted from them (labelled as full ). Results obtained using 84 regions are better than full atlases. This is probably due to the fact that small regions can induce spurious correlations in the connectivity matrices. Classification accuracy scores using overlapping atlases and their non overlapping counterpart. Point above the identity line support the idea that fuzzy overlapping maps are better at predicting behavioral variable than non-overlapping.

		All subjects	Biggest sites		Right-handed males	Right-hand males, 9-18	Right-hand males, 9-18, 3 sites
	Sites	ASD	TC	ASD	TC	ASD	TC	ASD	TC	ASD	TC
	CALTECH	5	10
	CMU	6	5
	KKI	12	21	12	21	7	15	6	12
	LEUVEN 1	14	14	14	14	13	13	1	1
	LEUVEN 2	12	16	12	16	7	12	7	12
	MAX MUN	19	27	19	27	16	25	5	3
	NYU	74	98	74	98	60	70	36	40	36	40
	OHSU	12	13			10	13	9	11
	OLIN	14	14	14	14	8	11		7
	Inter-site Figure A4: Comparison of prediction results obtained using feature selection or sparse methods. When selecting a reduced number of features, using either ANOVA to select 10% of the fea-tures, or sparsity inducing methods, we observe that the prediction accuracy drops compared to a classification using all features. This can be due to the fact that the connectivity matrices represent mixed 20 18 11 10 7 10 7 10 8 9 7 6 18 24 12 14 29 24 28 24 28 24 9 8 9 8 9 8 23 34 22 30 22 30 Figure A2: Intra-site PITT 24 26 24 26 SBL 12 14 SDSU 8 19 8 19 STANFORD 12 13 TRINITY 19 25 19 25 UCLA 1 37 27 37 27 UCLA 2 11 10 11 10 UM 1 34 52 34 52 50% 55% 60% 65% SVC-2 Ridge classifier ANOVA + SVC-2 ANOVA + Ridge classifier SVC-1 Lasso Without feature selection With feature selection UM 2 13 21 13 21 12 19 12 17 12 17 USM 43 24 43 24 37 22 12 7 YALE 22 19 22 19 8 10 7 10 Figure A3: 45% Total 403 468 356 413 292 337 198 222 107 119

features that should be handled together rather than seaprately. This may also be due to global effects in the connectivity: A global hypoconnectivity, as already observed in ASD patients, cannot be captured by 1 regularized classifiers.

Table A2 :

 A2 Count of subjects per site and per subsamples.

	Cross	Atlas estimator	Right handed	Right handed	Right handed	Biggest sites	All subjects
	validation		males 9-18 yo, 3	males 9-18 yo	males
			sites
		MSDL	66.6% ± 5.4%	65.8% ± 5.9%	65.7% ± 4.9%	67.9% ± 1.9%	66.9% ± 2.7%
		Yeo	60.9% ± 7.5%	62.3% ± 3.6%	64.7% ± 3.1%	64.7% ± 2.9%	66.9% ± 3.0%
	Intra-site	Harvard Oxford ICA	63.6% ± 2.7% 61.3% ± 8.7%	62.1% ± 4.6% 62.3% ± 3.0%	64.8% ± 5.1% 61.6% ± 4.3%	64.8% ± 2.5% 65.2% ± 3.1%	66.4% ± 2.7% 62.0% ± 3.4%
		K-Means	62.6% ± 5.9%	61.4% ± 5.4%	61.5% ± 3.0%	61.3% ± 3.4%	65.1% ± 1.8%
		Ward	63.2% ± 6.3%	60.1% ± 6.0%	62.2% ± 4.9%	64.6% ± 2.7%	63.7% ± 3.4%
		MSDL		68.3% ± 7.6%	63.4% ± 6.3%	68.7% ± 9.3%	66.8% ± 5.4%
		Yeo		69.7% ± 8.9%	64.5% ± 10.3%	61.4% ± 7.9%	61.3% ± 7.2%
	Inter-site	Harvard Oxford ICA		68.1% ± 9.0% 63.1% ± 9.9%	65.1% ± 5.8% 62.5% ± 7.8%	62.4% ± 5.4% 65.0% ± 4.8%	63.6% ± 6.2% 60.9% ± 5.2%
		K-Means		62.8% ± 13.9%	61.5% ± 8.1%	61.9% ± 10.1%	60.3% ± 4.8%
		Ward		62.4% ± 11.7%	59.8% ± 6.9%	63.4% ± 5.4%	63.1% ± 4.0%

Table A3 :

 A3 Average accuracy scores (along with their standard deviations) for top performing pipelines. This table summarizes the scores for pipelines with best prediction accuracy for each atlas and subset using intra-site or inter-site prediction. Best results are shown in bold.

	Cross	Atlas estimator	Right handed	Right handed	Right handed	Biggest sites	All subjects
	validation		males 9-18 yo, 3	males 9-18 yo	males
			sites
		MSDL	73.6% ± 11.7%	76.5% ± 25.1%	62.4% ± 21.3%	59.3% ± 8.0%	61.3% ± 4.8%
		Yeo	68.2% ± 36.5%	63.0% ± 20.2%	58.5% ± 4.7%	59.0% ± 5.1%	59.6% ± 5.6%
	Intra-site	Harvard Oxford ICA	67.7% ± 23.3% 69.1% ± 12.1%	61.4% ± 10.8% 65.7% ± 24.5%	59.5% ± 5.3% 64.4% ± 7.9%	57.4% ± 7.8% 57.3% ± 4.2%	59.5% ± 6.0% 61.3% ± 6.5%
		K-Means	69.1% ± 12.3%	68.9% ± 22.7%	61.4% ± 5.3%	53.3% ± 5.9%	60.4% ± 6.6%
		Ward	66.4% ± 20.9%	63.2% ± 21.4%	62.7% ± 5.8%	53.3% ± 4.9%	59.8% ± 5.9%
		MSDL		72.7% ± 19.0%	67.5% ± 10.8%	67.9% ± 9.0%	66.6% ± 14.7%
		Yeo		65.4% ± 13.7%	67.1% ± 14.7%	57.2% ± 9.4%	59.1% ± 8.8%
	Inter-site	Harvard Oxford ICA		71.9% ± 22.5% 69.9% ± 17.1%	66.7% ± 21.2% 65.8% ± 13.0%	60.7% ± 12.1% 63.4% ± 12.2%	64.2% ± 9.0% 68.9% ± 6.7%
		K-Means		69.7% ± 13.2%	67.7% ± 13.7%	60.3% ± 15.2%	67.7% ± 6.2%
		Ward		66.0% ± 12.1%	63.3% ± 16.9%	59.5% ± 9.3%	65.4% ± 9.7%

Table A4 :

 A4 Average specificity scores (along with their standard deviations) for top performing pipelines. This table summarizes the scores for pipelines with best prediction accuracy for each atlas and subset using intra-site or inter-site prediction. Best results are shown in bold.

	Cross	Atlas estimator	Right handed	Right handed	Right handed	Biggest sites	All subjects
	validation		males 9-18 yo, 3	males 9-18 yo	males
			sites
		MSDL	85.2% ± 14.6%	86.0% ± 12.9%	88.1% ± 12.6%	84.9% ± 11.7%	90.9% ± 7.5%
		Yeo	82.0% ± 13.2%	82.1% ± 8.4%	84.3% ± 12.6%	81.1% ± 6.1%	79.7% ± 3.4%
	Intra-site	Harvard Oxford ICA	83.6% ± 14.3% 80.4% ± 16.4%	79.1% ± 20.4% 77.7% ± 11.8%	85.9% ± 10.6% 88.4% ± 8.2%	84.4% ± 3.0% 81.5% ± 8.2%	82.3% ± 2.7% 85.1% ± 4.4%
		K-Means	74.0% ± 13.1%	78.8% ± 20.3%	79.4% ± 19.0%	83.2% ± 5.7%	80.0% ± 12.2%
		Ward	78.0% ± 11.0%	77.0% ± 11.3%	80.0% ± 16.2%	80.4% ± 3.6%	80.4% ± 4.3%
		MSDL		84.2% ± 10.5%	85.4% ± 12.5%	85.5% ± 13.6%	100.0% ± 0.0%
		Yeo		80.4% ± 12.3%	81.5% ± 16.0%	76.5% ± 16.3%	80.0% ± 15.1%
	Inter-site	Harvard Oxford ICA		69.7% ± 16.5% 71.2% ± 16.2%	71.9% ± 10.9% 72.5% ± 7.6%	73.3% ± 13.6% 75.5% ± 12.9%	88.7% ± 18.9% 93.1% ± 15.0%
		K-Means		75.4% ± 18.3%	67.2% ± 14.3%	68.4% ± 16.8%	94.3% ± 18.1%
		Ward		70.3% ± 17.0%	63.2% ± 11.8%	74.1% ± 10.8%	82.6% ± 22.9%
	Table

See http://fcon_1000.projects.nitrc.org/indi/abide/ for specific information

http://preprocessed-connectomes-project.github.io/abide

For the sake of simplicity, we report only the 3 main methods here. However, we also investigated random forests and Gaussian naive Bayes, that gave low prediction performance. Feature selection also led to worse results, both with univariate ANOVA screening and with sparse models such as Lasso or SVC-1 -see FigureA4for more details. Logistic regression gave results similar to SVC.

Software versions: python 2.7, scikit-learn 0.17.1, scipy 0.14.0, numpy 1.9.0, nilearn 0.1.5.

Table A6 lists subjects of each subsample. The lists of subjects for each cross-validation fold are available at https://team.inria. fr/parietal/files/2016/04/cv_abide.zip

Note that the criterion used to select the data-point used for the post-hoc analysis, i. e. the 10% best performing, is independent from the factors studied in this analysis, as we select the 10% best performing for each value of these factors. Hence, the danger of inflating effects by a selection[START_REF] Vul | Puzzlingly high correlations in fmri studies of emotion, personality, and social cognition[END_REF] does not apply to the post hoc analysis.

-regularized classifiers, which induce sparsity. In addi-

tion 2 -penalization is rotationally-invariant, which means that it is not sensitive to an arbitrary mixing of features. Such mixing may happen if a functional brain module sits astride several regions. Thus, a possible explanation for the good performance of 2 -penalization is that it does not need an optimal choice of regions, perfectly aligned with the underlying functional modules.

Note that we experimented with a variety of more sophisticated covariance estimators, including GraphLasso, though without improved results.

Acknowledgments. We acknowledge funding from the NiConnect project and the SUBSample project from the DIGITEO Institute, France. The effort of Adriana Di Martino was partially supported by the NIMH grant 1R21MH107045-01A1. Finally, we thank the anonymous reviewers for their feedback, as well as all sites and investigators who have worked to share their data through ABIDE.

Appendix A. Time series extraction

In order to extract time series from overlapping brain maps, we use an ordinary least square approach (OLS). Let Y ∈ R n×p be the original subject signals of p voxels over n temporal scans and V ∈ R k×p the atlas of k maps containing p voxels. We estimate U ∈ R n×k , the set of time series proper to each region, using an OLS:

In order to remove the confounds, we orthogonalize our signals with regards to them. We do this by removing the the signals projected on an orthonormal basis of the confounds. Let C ∈ R n×c be a matrix of c confounds:

Covariance matrix estimation

The empirical covariance matrix is not a good estimator when n < p. In particular, it can leads to the estimation of matrices that cannot be inverted. This problem can be avoided by using shrinkage, a method to pull the extreme values of the covariance matrix toward the central ones. It is done by finding the 2 -penalized maximum likelihood estimator of the covariance matrix:

We use the Ledoit-Wolf approach [START_REF] Ledoit | A well-conditioned estimator for largedimensional covariance matrices[END_REF] to set parameter α to minimize the Mean Square Error between the real and the estimated covariance matrix.

Appendix C. Classifiers

Linear classifiers try to find an hyperplane splitting the sample in the input space in order to classify them into two categories.

Support Vector Classification SVC

SVC aims at maximizing the margin, i. e. the distance to the closest training sample and the hyperplane. We use the classical squared hinge loss and thus minimize: ). Each bar represents the impact of the corresponding option on the prediction accuracy, relatively to the mean prediction. This effect is measured via a full-factorial ANOVA, analyzing the contribution of each step in a linear model. Each step of the pipeline is considered as a categorical variable. Error bars give the 95% confidence interval. MSDL atlas extraction method gives significantly better results while reference atlases are slightly better than the mean. Among all matrix types, tangent embedding is the best on all ABIDE subsets. Finally, l 2 regularized classifiers are better than l 1 regularized.