Kamalaker Dadi

Mehdi Rahim

Alexandre Abraham

Darya Chyzhyk

Michael Milham

Bertrand Thirion

Gaël Varoquaux

Benchmarking functional connectome-based predictive models for resting-state fMRI

Keywords: Resting-state fMRI, Functional connectomes, Predictive modeling, Classification, Population study

come

Introduction

Resting-state functional Magnetic Resonance Imaging (rest-fMRI), based on the analysis of brain activity without specific task, has become a tool of choice to probe human brain function in healthy and diseased populations. As it can easily be acquired in many different individuals, rest-fMRI is a promising candidate for markers of brain function [START_REF] Biswal | Toward discovery science of human brain function[END_REF][START_REF] Greicius | Resting-state functional connectivity in neuropsychiatric disorders[END_REF]. This has lead to the rise of large-scale rest-fMRI data collections, such as the human connectome project [START_REF] Van Essen | The wu-minn human connectome project: an overview[END_REF] or ABIDE [START_REF] Martino | The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism[END_REF]. Larger datasets bring increased statistical power [START_REF] Elliott | The UK biobank sample handling and storage protocol for the collection, processing and archiving of human blood and urine[END_REF], and many population-imaging studies use rest-fMRI to relate brain imaging to neuropathologies or other behavior and population phenotypes [START_REF] Miller | Multimodal population brain imaging in the UK biobank prospective epidemiological study[END_REF][START_REF] Dubois | Building a science of individual differences from fmri[END_REF]. These efforts build biomarkers from rest-fMRI with predictive models [START_REF] Woo | Building better biomarkers: brain models in translational neuroimaging[END_REF]. * Corresponding author A functional connectome -characterizing the network structure of the brain [START_REF] Sporns | The human connectome: a structural description of the human brain[END_REF]-can be extracted from functional interactions in rest-fMRI data [START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF]. The weights of the corresponding brain functional connectome are used to characterize individual subjects behavior, cognition, and mental health [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] Milazzo | Identification of moodrelevant brain connections using a continuous, subject-driven rumination paradigm[END_REF][START_REF] Smith | A positive-negative mode of population covariation links brain connectivity, demographics and behavior[END_REF][START_REF] Miller | Multimodal population brain imaging in the UK biobank prospective epidemiological study[END_REF][START_REF] Colclough | The heritability of multi-modal connectivity in human brain activity[END_REF][START_REF] Dubois | A distributed brain network predicts general intelligence from resting-state human neuroimaging data[END_REF], aging [START_REF] Liem | Predicting brain-age from multimodal imaging data captures cognitive impairment[END_REF] as well as brain pathologies [START_REF] Drysdale | Resting-state connectivity biomarkers define neurophysiological subtypes of depression[END_REF][START_REF] Abraham | Deriving reproducible biomarkers from multi-site resting-state data: An autism-based example[END_REF][START_REF] Ng | Distinct alterations in parkinson's medicationstate and disease-state connectivity[END_REF].

Machine-learning pipelines are key to turning functional connectomes into biomarkers that predict the phenotype of interest [START_REF] Woo | Building better biomarkers: brain models in translational neuroimaging[END_REF]. On rest-fMRI, such a pipeline typically comprises of 3 crucial steps as depicted in Figure 1, linking functional connectomes to the target phenotype [START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF][START_REF] Craddock | Connectomics and new approaches for analyzing human brain functional connectivity[END_REF]. Yet, there exist many variations of this prototypical pipeline, even for classification from edge-weights of brain functional connectomes, as revealed by reviews of the field

RS-fMRI

Diagnosis

Connectivity

Parameterization

Supervised Learning Defining Brain ROIs 1 2 3 [START_REF] Wolfers | From estimating activation locality to predicting disorder: a review of pattern recognition for neuroimaging-based psychiatric diagnostics[END_REF][START_REF] Arbabshirani | Single subject prediction of brain disorders in neuroimaging: Promises and pitfalls[END_REF][START_REF] Brown | Machine learning on human connectome data from mri[END_REF]. These various choices have a sizable impact on the accuracy of population studies, and are seldom discussed [START_REF] Carp | On the plurality of (methodological) worlds: Estimating the analytic flexibility of fMRI experiments[END_REF]. The cost of such analytical variation is twofold. First, it puts the burden on the practitioner to explore many options and make choices without systematic guidance. Second, methods variations create researchers degrees of freedom [START_REF] Simmons | False-positive psychology: Undisclosed flexibility in data collection and analysis allows presenting anything as significant[END_REF] that can compromise the measure of the prediction accuracy of biomarkers [START_REF] Varoquaux | Cross-validation failure: small sample sizes lead to large error bars[END_REF]. Guidelines on optimal modeling choices are thus of great value for rest-fMRI biomarker research.

Here, we perform a systematic benchmark of common choices for the different steps of the functional connectomebased classification pipeline. To outline the preferable strategies, we analyze the prediction accuracy across 6 different cohorts, with different clinical questions and one psychological trait, different sample sizes, and prediction problems of different difficulties. While best model choice may vary depending on the prediction task, our benchmarks outline some trends. Specifically, we explore the following analytical choices:

• How should nodes be chosen: via pre-defined atlases, or data-driven approaches? How many nodes are needed for brain-imaging based diagnosis? Should nodes be distributed brain networks or regions of interest (ROIs)? • How should weights of brain functional connectomes be represented: via correlations, partial correlation, or more complex models capturing the geometry of covariance matrices? • What classifiers should be used for machine learning on weights of brain functional connectomes? Should linear or non-linear models be preferred? Should sparse or nonsparse models be used? With or without feature selection? Besides these main questions, we did additional experiments on preprocessing strategies -studying the effect of band-pass filtering and global signal regression-and on covariance estimators, by comparing sparsity-inducing to classical shrinkage.

The paper is organized as follows: we first review current practices and methods used to-date for prediction of psychiatric diseases from weights of brain functional connectomes. Then, we present the different choices that we benchmark for the steps of classification pipelines and describe these methods. Finally, we report our experimental results and the trends that they reveal.

Methods: functional connectome-classification pipeline

Figure 1 shows the standard rest-fMRI classification pipeline that we consider.

A brief review of current practices: functional connectome-based predictive methods

We first survey methods used for prediction studies based on three extensive reviews: [START_REF] Wolfers | From estimating activation locality to predicting disorder: a review of pattern recognition for neuroimaging-based psychiatric diagnostics[END_REF]; [START_REF] Arbabshirani | Single subject prediction of brain disorders in neuroimaging: Promises and pitfalls[END_REF]; [START_REF] Brown | Machine learning on human connectome data from mri[END_REF]. From these reviews, 27 studies used rest-fMRI and gave good classification scores. Below, we briefly outline the choices in the different pipeline step used (see Table A2 in the appendix for the full list).

Definition of brain ROIs. Studies define ROIs to extract signals with a variety of approaches:

• balls1 of radius varying from 5mm to 10mm centered at coordinates from the literature [START_REF] Dosenbach | Prediction of individual brain maturity using fMRI[END_REF][START_REF] Power | Functional network organization of the human brain[END_REF]; • reference anatomical atlases such as AAL [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], sulci-based atlases [START_REF] Perrot | Joint Bayesian Cortical Sulci Recognition and Spatial Normalization[END_REF][START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on mri scans into gyral based regions of interest[END_REF], or connectivity-based cortical landmarks [START_REF] Zhu | Connectome-scale assessments of structural and functional connectivity in MCI[END_REF]; • data-driven approaches based on k-means or Ward clustering, as well as Independent Component Analysis (ICA) approaches [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] or dictionary learning [START_REF] Abraham | Extracting brain regions from rest fMRI with totalvariation constrained dictionary learning[END_REF]. The number of nodes used was typically around 100, but ranged from dozens to several hundreds.

Representation of brain functional connectomes. Studies define functional interactions from second-order statistics -based on signal covariance-using Pearson's correlation or partial correlations estimated mostly either with the maximum-likelihood formula for the covariance or the Ledoit-Wolf shrinkage covariance estimator [START_REF] Ledoit | A well-conditioned estimator for largedimensional covariance matrices[END_REF][START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF][START_REF] Brier | Partial covariance based functional connectivity computation using ledoit-wolf covariance regularization[END_REF]. Partial correlation between nodes is useful to rule out indirect effects in the correlation structure, but calls for shrunk estimates [START_REF] Smith | Network modelling methods for fMRI[END_REF][START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF]. Mathematical arguments have also led to representations tailored to the manifold-structure of covariance matrices (Varoquaux et al., 2010a;[START_REF] Ng | Transport on Riemannian Manifold for Functional Connectivity-based Classification[END_REF][START_REF] Dodero | Kernel-based classification for brain connectivity graphs on the riemannian manifold of positive definite matrices[END_REF][START_REF] Colclough | The heritability of multi-modal connectivity in human brain activity[END_REF]. We benchmark the simplest of these, a tangent representation of the manifold which underlies the more complex developments (see Appendix A for a quick introduction to this formalism).

Classifiers used for prediction. Many different classifiers have been used, whether linear or non-linear, sparse or non-sparse, optionally with prior feature selection. See Table A2 for the comprehensive list of classifiers used in these studies.

Finally, beyond the prototypical pipeline exposed above, some studies employ complex-graph network modeling approaches -e.g. network modularity or centrality [START_REF] Rubinov | Weight-conserving characterization of complex functional brain networks[END_REF])- [START_REF] Wolfers | From estimating activation locality to predicting disorder: a review of pattern recognition for neuroimaging-based psychiatric diagnostics[END_REF][START_REF] Arbabshirani | Single subject prediction of brain disorders in neuroimaging: Promises and pitfalls[END_REF][START_REF] Brown | Machine learning on human connectome data from mri[END_REF] These approaches are seldom combined with supervised learning. Indeed, graph-theory metrics capture well global aspects of brain connectivity, but do not lend themselves well to tuning to connections in specific subnetworks [START_REF] Hallquist | Graph theory approaches to functional network organization in brain disorders: A critique for a brave new small-world[END_REF]. Here, we focus on machine-learning methods that extract discriminant connections; as such we do not study graph-theoretical approaches.

The current practice is very diverse, without standard modeling choices. To open the way toward informed decisions, we explore popular variants of the classic machinelearning pipeline to predict on connectomes. We measure the impact of choices at each step on prediction for diverse targets across multiple datasets. We detail below the specific modeling choices included in our benchmarks.

Definition of brain regions of interest (ROIs)

For functional connectomes, the hypothesis is that the definition of ROIs should capture well the relevant functional units [START_REF] Smith | Network modelling methods for fMRI[END_REF]. We study both anatomically and functionally defined reference brain atlases, as well as data-driven methods that define ROIs from the data at hand. ROI selection is a difficult choice, as the optimal may vary for different conditions or pathologies.

A selection of pre-defined atlases. We consider four standard atlases, of which two are structural atlases: i) Automated Anatomical Labeling (AAL) [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], a structural atlas with 116 ROIs defined from the anatomy of a reference subject, ii) 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 probabilistic atlas of anatomical structures, contains of 48 cortical & 11 subcortical ROIs in each hemisphere, ie 118 ROIs in total. We also include two functional atlases: iii) Bootstrap Analysis of Stable Clusters (BASC) [START_REF] Bellec | Multi-level bootstrap analysis of stable clusters in resting-state fMRI[END_REF], a multi-scale functional atlas built with clustering on rest-fMRI, coming with different {36, 64, 122, 197, 325, 444} numbers of ROIs; iv) Power, a coordinate-based atlas consisting of 264 coordinates used to position balls of 5mm radius [START_REF] Power | Functional network organization of the human brain[END_REF]. For an additional set of benchmarks, on larger data, we use only pre-computed regions. For a pre-computed functional atlas with dictionary learning, we use an atlas2 computed by Mensch et al. (2016a) with a very scalable sparse dictionary-learning algorithm on the HCP900 dataset [START_REF] Van Essen | The human connectome project: A data acquisition perspective[END_REF]. This algorithm, MODL (massive online dictionary learning), solves the 1 dictionary-learning problem with an algorithm fast on very large datasets that converges to the same solution as standard on-line solvers [START_REF] Mensch | Stochastic subsampling for factorizing huge matrices[END_REF].

A selection of data-driven methods. We consider four popular data-driven methods to extract brain ROIs from intrinsic brain activity [START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF][START_REF] Kahnt | Connectivity-based parcellation of the human orbitofrontal cortex[END_REF][START_REF] Thirion | Which fMRI clustering gives good brain parcellations? Name[END_REF][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][START_REF] Abraham | Extracting brain regions from rest fMRI with totalvariation constrained dictionary learning[END_REF]. We choose to define ROIs using two clustering methods: i) K-Means [START_REF] Hastie | The elements of statistical learning[END_REF], and ii) hierarchical agglomerative clustering using Wards algorithm [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF] with spatial connectivity constraints [START_REF] Michel | A supervised clustering approach for fMRIbased inference of brain states[END_REF]; and two linear decomposition methods: iii) Canonical Independent Component Analysis (GroupICA or CanICA) (Varoquaux et al., 2010c), iv) Dictionary Learning -1 (DictLearn) (Mensch et al., 2016b).

Dimension selection in data-driven atlases. For clustering methods, we extract brain atlases with a varying number of ROIs in dim = {40, 60, 80, 100, 120, 150, 200, 300}. With linear decomposition methods i.e. CanICA and DictLearn, we explore the following number of components: dim = {40, 60, 80, 100, 120}3 .

For each data-driven method, we learn brain ROIs on the training set only, to avoid possible overfit [START_REF] Abraham | Deriving reproducible biomarkers from multi-site resting-state data: An autism-based example[END_REF]. In a cross-validation loop, for each split, we define the brain ROIs on a training set and use the atlases to learn connectivity patterns for prediction. We also applied additional Gaussian smoothing of 6mm on preprocessed rest-fMRI datasets for all data-driven methods prior to learning brain ROIs to enhance the region extraction step.

Nodes formed local regions or distributed networks?. Current practices in functional connectomics includes defining nodes as local regions of the brain [START_REF] Shirer | Decoding subject-driven cognitive states with whole-brain connectivity patterns[END_REF]; Craddock et al., 2012), or as full distributed functional networks that may include several regions [START_REF] Smith | A positive-negative mode of population covariation links brain connectivity, demographics and behavior[END_REF][START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF]. We consider both approaches: using the distributed networks, or breaking them up in regions with a segmentation step to separate out regions (Abraham et al., 2014a). For example, a bi-hemispheric brain networks is separated into regions, one in each hemisphere.

We use a Random-Walker based extraction of regions from the brain networks obtained by CanICA and DictLearn as proposed in Abraham et al. (2014a). By contrast, for K-Means and BASC, we simply break out clusters in their connected components. During this procedure, we remove spurious regions of size < 1500mm 3 . Figure 2 shows an example of the set of brain regions obtained from the various data-driven methods on the ADNI rest-fMRI data.

Connectivity parametrization

We extract representative time series for each node. For signal extraction, we explored several denoising strategies to account for non-neural artifacts: with or without low-pass filtering or global signal mean regression (details in Appendix B). To estimate functional connectomes efficiently, we use the Ledoit-Wolf regularized shrinkage estimator [START_REF] Ledoit | A well-conditioned estimator for largedimensional covariance matrices[END_REF][START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF][START_REF] Brier | Partial covariance based functional connectivity computation using ledoit-wolf covariance regularization[END_REF], which gives a closed form expression for the shrinkage parameter. This estimator yields well-conditioned estimators despite the variation in length of time series across rest-fMRI datasets. We also explored non-regularized and sparse estimator for the covariance (see Appendix H.2). With this covariance structure, we study three different parametrizations of functional interactions: full correlation, partial correlation [START_REF] Smith | Network modelling methods for fMRI[END_REF][START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF] and the tangent space of covariance matrices. The latter is less frequently used but has solid mathematical foundations and a variety of groups have reported good decoding per-formances with this framework (Varoquaux et al., 2010a;[START_REF] Barachant | Classification of covariance matrices using a riemannian-based kernel for bci applications[END_REF][START_REF] Ng | Transport on Riemannian Manifold for Functional Connectivity-based Classification[END_REF][START_REF] Dodero | Kernel-based classification for brain connectivity graphs on the riemannian manifold of positive definite matrices[END_REF][START_REF] Qiu | Manifold learning on brain functional networks in aging[END_REF][START_REF] Rahim | Population-shrinkage of covariance to estimate better brain functional connectivity[END_REF][START_REF] Wong | Riemannian regression and classification models of brain networks applied to autism[END_REF]. We compared two variants, using as a reference point the Euclidean mean (Varoquaux et al., 2010a) or the geometric mean [START_REF] Ng | Transport on Riemannian Manifold for Functional Connectivity-based Classification[END_REF]; in both cases we rely on Nilearn implementation (Abraham et al., 2014b). Note that computing partial correlation or tangent space require inverting covariance matrices, hence these must be well conditioned. Non regularized covariance estimation is thus not useable for these parametrizations.

For each parametrization, we vectorize the functional connectome, using the lower triangular part of the connectomes matrix for classification. Additionally, on the ACPI dataset, we considered the ADHD status of the subjects as a variable of non interest and regressed it out in this second-level analysis, as we were interested in predicting the consumption of Marijuana.

Supervised learning: Classifiers

The final step of our pipeline predicts a binary phenotypic status from connectivity features extracted from previous step. We consider several linear and non-linear classifiers for prediction i.e. both sparse and non-sparse methods. For non-linear methods, we consider Nearest Neighbors (K-NN) [START_REF] Cover | Nearest neighbor pattern classification[END_REF] with K=1 and Euclidean distance metric, Gaussian Naïve Bayes (GNB) and Random Forests Classifier (RF) [START_REF] Breiman | Random forests[END_REF]. For linear classifiers we consider sparse 1 regularization4 for Support Vector Classification (SVC), and Logistic Regression [START_REF] Hastie | The elements of statistical learning[END_REF]. For non-sparse linear classifiers -i.e. 2 regularization-we consider Ridge classification, SVC, Logistic regression. For SVC, we also considered 10% feature screening with univariate ANOVA. With regards to the regularization parameter (eg soft margin parameter in SVC), we use the default C = 1 or α = 1, which has been found to be a good default [START_REF] Varoquaux | Assessing and tuning brain decoders: cross-validation, caveats, and guidelines[END_REF].

Experimental study

To benchmark the various predictive-modeling choices, we apply the functional connectome-classification pipeline on five publicly-available rest-fMRI datasets. We study prediction from functional connectomes of various clinical outcomes -neuro-degenerative and neuro-psychiatric disorders, drug abuse impact, fluid intelligence. We focus on binary classification problems, predicting a phenotypic target between two groups. We use the following datasets, summarized in Table 1: 1. COBRE, Center for Biomedical Research Excellence5 , comprising rest-fMRI data to study schizophrenia and bipolar disorder [START_REF] Calhoun | Exploring the psychosis functional connectome: Aberrant intrinsic networks in schizophrenia and bipolar disorder[END_REF]. We focus on predicting schizophrenia diagnosis versus normal control.

2. ADNI, the Alzheimer's Disease Neuroimaging Initiative 6 database studies neuro-degenerative diseases [START_REF] Mueller | The alzheimers disease neuroimaging initiative[END_REF]. We focus on using rest-fMRI to discriminate individuals with Mild Cognitive Impairment (MCI) from individuals diagnosed with Alzheimer's Disease (AD).

3. ADNIDOD, funded by the US Department of Defense (DoD) to study brain aging in Vietnam War Veterans 7 , includes rest-fMRI data of individuals with post-traumatic stress disorders (PTSD) or brain traumatic injuries. We focus on discriminating PTSD condition from normal controls.

4. ACPI, Addiction Connectome Preprocessed Initiative 8 , a longitudinal study to investigate the effect of cannabis use among adults with a childhood diagnosis of ADHD. In particular we use readily-preprocessed rest-fMRI data from Multimodal treatment study of Attention Deficit Hyperactivity Disorder (MTA). We attempt to discriminate whether individuals have consumed marijuana or not.

5. ABIDE, Autism Brain Imaging Data Exchange database investigates the neural basis of autism (Di [START_REF] Martino | The autism brain imaging data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism[END_REF]. We use the data from Preprocessed Connectome Project [START_REF] Craddock | The neuro bureau preprocessing initiative: open sharing of preprocessed neuroimaging data and derivatives ¡br /¿[END_REF] to discriminate individuals from Autism Spectrum Disorder from normal controls.

6. HCP, 9 Human Connectome Project contains imaging and behavioral data of healthy subjects (Van Essen et al., 2013). We use preprocessed rest-fMRI data from HCP900 release [START_REF] Van Essen | The human connectome project: A data acquisition perspective[END_REF] to discriminate individuals from high IQ and low IQ. We used HCP rs-fMRI datasets to probe a different setting: data with longer acquisitions. Due to the data size, we limit the benchmarks here to pre-computed atlases.

3.1. rest-fMRI data processing: softwares and related Data preprocessing. We preprocess COBRE, ADNI, and ADNIDOD. We use a standard protocol that includes: motion correction, fMRI co-registration to T1-MRI, normalization to the MNI template using SPM12 Exclusion criteria. We not only exclude subjects based on visual inspection of preprocessed data, but also subjects that do not fall into binary classification groups, eg we removed subjects who had both bipolar disorder and schizoaffective groups from COBRE samples. For HCP, we select the subjects with single session and phase encoding in a left-to-right (LR) direction. Out of these selected subjects, we discriminate the low IQ from the high IQ individuals, where the data are split in 3 according to quantiles 0.333 and 0.666, and the subjects in the middle group are excluded to make the prediction easier in a binary classification setup (see Table 1 for numbers of subjects included in the analysis).

Cross validation and error measure. We perform crossvalidation (CV) by randomly shuffling and splitting each dataset over 100 folds, forming two sets of subjects: 75% for training the classifier and learning brain atlases with data-driven models and the remaining 25% for testing on unseen data [START_REF] Varoquaux | Assessing and tuning brain decoders: cross-validation, caveats, and guidelines[END_REF]. We create stratified folds, preserving the ratio of samples between groups. For each split, we measure the Area Under the Curve (AUC) from the Receiver Operating Characteristics (ROC) curve: 1 is a perfect prediction and .5 is chance. The final prediction scores in AUC (> 120k scores) are used to measure the impact of various choices in our prediction pipeline outlined below in results section.

Computations and implementation. Our experimental study consists of more than 240 types of pipelines (8 at-lases × 3 connectivity measures × 10 classifiers, plus some variants such as 3 filtering options and 3 covariance estimator options). These pipelines were run on each of 5 datasets for 100 CV folds. As a result, there are more than 500 000 pipeline fits, from the raw data to the supervised step, a heavy computational load. Technically, we rely on efficient implementations open-source scientific computing packages using Python 2.7: Nilearn v0.3 (Abraham et al., 2014b) to define brain atlases, extract representative timeseries and timeseries confounds regression, and build connectivity measures. All machine-learning methods used for prediction i.e., classifiers and cross-validation are implemented with scikit-learn v0.18.1 [START_REF] Pedregosa | Scikitlearn: Machine learning in Python[END_REF].

For visualization, we rely on Nilearn for brain-related figures while matplotlib is used [START_REF] Hunter | Matplotlib: A 2d graphics environment[END_REF] for generating other figures.

Results: benchmarks of pipeline choices

We now outline which modeling choices have an important impact on predicting over diverse phenotypes from all rest-fMRI datasets.

We report in Table 2 the AUC scores obtained for all rest-fMRI datasets. The scores reported in the table are simplified to the optimal choice selection at each step in the pipeline which showed significant impact. These optimal choice of steps are discussed in following sections.

Impact of methodological choices. We study the prediction score of each pipeline relative to the mean across pipelines on each fold. This relative measure discards the variance in scores due to folds or datasets. From these relative prediction scores, we study the impact of the choice of each step in the prediction pipeline: choice of classifiers, connectivity parametrizations, and brain atlases. This is a multifactorial set of choices and there are two points of view on the impact of a choice for a given step. First, the impact of the choice for one step may be considered when the other steps are optimal, or close to optimal. Second, the impact of one step may be considered for all other choices for the other steps -marginally on the choice of other steps. Empirically, the two scenario lead to similar conclusions. In the following figures, we study the first Accuracy scores reported correspond to optimal choices in functional connectivity prediction pipeline: brain regions defined with regions using DictLearn, connectivity matrices parametrized by their tangent-space representation, and an 2 -regularized logistic regression as a classifier, as discussed below.

Best prediction is achieved with schizophrenia vs control discrimination task on COBRE dataset at 86.2% (median).

Figure 3: Impact of classifier choices on prediction accuracy, for all rest-fMRI datasets and all folds. For each classifier choice, only the top third highest performing scores are represented when varying the modeling choices for other steps in the pipeline: brainregion definition and connectivity parametrization. Figure A2 gives all the data points, not limited to good choices in the overall pipeline.

Overall, 2 -regularized linear classifiers perform better, with a slight lead for 2 logistic regression. The box plot gives the distribution across folds (n=100) and datasets (denoted by markers) of prediction score for a given choice (classifier) relative to the mean across all choices (regions-definition and connectivity parametrizations, classifiers). The box displays the median and quartiles, while the whiskers give the 5 th and 95 th percentiles.

situation, focusing on "good choices": given a choice for one step, we report data for top third highest performing scores (quantiles 0.666) for the choices in the other steps. Appendix C gives results for all scores, hence studying one choice, marginally upon the others.

Choice of classifier

Figure 3 summarizes the performances of classifiers on prediction scores for all rest-fMRI datasets. The results display a certain amount of variance across folds and datasets (i.e., prediction targets). However, they show that non-sparse ( 2 -regularized) linear classifiers perform better, with a slight lead for logistic-2 . Using non-linear classifiers does not appear useful; neither does sparsity. The results in Figure 3

Choice of connectivity parameterization

Figure 4 summarizes the impact of covariance matrix parametrization on the relative prediction scores for all rest-fMRI datasets. Tangent-space parametrization tends Figure 4: Impact of connectivity parameterization on prediction accuracy, for all rest-fMRI datasets and folds. For each parametrization choice, only the top third highest performing scores are represented when varying the modeling choices for other steps in the pipeline: brain-region definition and classifier. Figure A3 gives all the data points, not limited to good choices in the overall pipeline. Prediction using tangent space based connectivity parameterization displays higher accuracy with relatively lower variance than using full or partial correlation. The box displays the median and quartiles, while the whiskers give the 5 th and 95 th percentiles.

to outperform full correlations or partial correlations. Indeed, it performs better on average, but also has less variance across datasets (prediction targets) or folds. Results are similar for simpler variant of the tangent-space parametrization relying on a simple Euclidean mean rather than the full geometric (Riemannian) -see Appendix A for more details. While scores in Figure 4 are conditional on a good choice for other pipeline parameters, Figure A3 gives results marginal to all choices. In both settings, connectivity matrices built with tangent space parametrization give an improvement compared to full or partial correlations.

Choice of regions definition method

To find the preferred approaches to define brain regions, we proceed in two steps. First, for each method, we find the dimensionality that gives the best prediction. This holds for the BASC atlas, that comes in various dimensionalities, and for data-driven region-definition methods, for which we vary the dimensionality. Second, we study the prediction accuracy for each approach at the optimal dimensionality. Best approach. Figure 5 summarizes the relative prediction performance of all choices of region-definition methods. While the systematic effects are small compared to the variance over the folds and the datasets, the general trend is that regions defined from functional data lead to better prediction than regions defined from anatomy. Using 1 dictionary learning to define regions from rest-fMRI data appears to be the best method, closely followed by ICA, which is also based on a linear decomposition model. Interestingly, BASC, an atlas pre-defined on unrelated rest-fMRI datasets using data-driven clustering Figure 5: Impact of region-definition method on prediction accuracy, for all rest-fMRI datasets and folds. For each regiondefinition choice, only the top third highest performing scores are represented when varying the modeling choices for other steps in the pipeline: classifier and connectivity parametrization. Figure A4 gives all the data points, not limited to good choices in the overall pipeline. Learning atlases from rest-fMRI data tends the prediction for all tasks. By contrast anatomical atlases perform poorly over diverse tasks. The box displays the median and quartiles, while the whiskers give the 5 th and 95 th percentiles.

technique, performs almost as well as the best regionsextraction method applied to the rest-fMRI data of interest. Unlike other pre-defined atlases, like Harvard Oxford or AAL, that lack some crucial functional regions. The BASC atlas [START_REF] Bellec | Multi-level bootstrap analysis of stable clusters in resting-state fMRI[END_REF] is readily available online, and is thus easy to apply to data. Figure 5 shows the impact of region-definition approach conditional on good choices in the other steps of the pipeline, however studying the impact of region-definition independently of other choices (Figure A4). Both comparisons highlight that defining regions from functional data gives the bestperforming pipelines, and that linear-decomposition methods are to be preferred.

Optimal dimensionality. The choice of the best dimensionality for each approach paints a less clear picture (Figure 6): a range of dimensionalities lead for good prediction for each method12 . We find that there is a very soft optimum: prediction reaches a plateau as the number of extracted networks increases, and then slowly decreases for some methods. To favor the most parsimonious model, in this paper we choose to work at the lower end of the plateau (red arrow on Figure 6): simpler models for better stability and statistical control. While this choice is not clear cut, the curves also suggest that, in a reasonable range, it does not have a large impact on prediction accuracy. Note that the dimensionality here corresponds to the number of networks, these are then broken up into separate regions. We find that the typical number of regions at the optimal is around 150 (Appendix E).

Localized regions or distributed networks. Nodes of the functional connectomes may be defined from localized regions, or the distributed networks that naturally arise from approaches such as ICA or dictionary learning. The choice of one over the other has little impact over prediction, though there is slight, non significant, benefit to using regions (Figure A5).

Larger datasets and pre-computed atlases

To investigate the consistency of analytics choices for higher-quality datasets, we perform extra benchmarks including the HCP data. As this data comprises much longer time-series, we restrict our analysis to pre-computed atlases, that alleviate computational costs. We share the resulting time-series and scripts to reproduce our analysis 13 .

Figure 7 summarizes the impact of method choice on the prediction accuracy for all six different cohorts. This experiment outline similar tradeoffs as the others: functional atlas pre-computed with dictionary learning (here MODL, from Mensch et al. (2016a)), tangent-space parametrization, and 2 -regularized classifiers are preferable. This experiment is not as systematic as the other, as a very large dataset like HCP would require much more 13 github.com/KamalakerDadi/benchmark rsfMRI prediction computing power to study region extraction14 . Yet, even for region-definition methods, it outlines similar trends than when tuning the regions to the data at hand.

Filtering, global signal, and covariance estimation

When extracting functional signal for connectivity modeling, there are many options to reject confounds, including temporal filtering or global signal mean regression [START_REF] Fox | The global signal and observed anticorrelated resting state brain networks[END_REF][START_REF] Murphy | The impact of global signal regression on resting state correlations: are anti-correlated networks introduced[END_REF][START_REF] Power | Spurious but systematic correlations in functional connectivity mri networks arise from subject motion[END_REF]. Also, to extract a structure reflecting well brain connectivity, it has been show that careful covariance estimation is useful and that sparse inverse covariance methods perform well [START_REF] Smith | Network modelling methods for fMRI[END_REF][START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF]. For both of these steps, our experiments for predictive modeling applications do not reveal clear preferences.

Different time-series filtering approaches (band-pass or global-signal regression) make no visible differences on prediction accuracy (Figure A9). A likely reason is that the supervised step can learn a predictor that is independent of the corresponding noise in the signal.

With regards to covariance estimation, we also investigate the empirical covariance (maximum likelihood estimator) and sparse inverse covariance (Appendix H). The empirical covariance can only be used to compute correlations -as partial correlation or tangent parametrization require an invertible covariance matrix-in which case it performs similarly as the Ledoit-Wolf estimator (see Figure A11). Sparse inverse covariance performs as well or worse than the Ledoit-Wolf estimator. This latter estimator is easier to use, as it is faster and does not require setting a regularization parameter. Learning discriminant connectivity patterns across conditions does not seem to require the same regularization -sparsity-as identifying the brain connectivity structure.

- MODL is a parcellation built using a form of Online dictionary learning. Restricting to pre-computed regions and adding a different dataset (HCP) gives results consistent with Figure 5, 4, and 3: best choices are regions defined functionally, with decomposition methods (MODL) followed by clustering methods (BASC), tangent-space parametrization of connectivity, and 2 -regularized logistic regression. The box displays the median and quartiles, while the whiskers give the 5 th and 95 th percentiles. Table A1 reports the corresponding absolute scores.

Discussion

An increasing amount of studies use predictive models on functional connectomes, for instance in populationimaging settings to relate brain activity to psychological traits or to build biomarkers of pathologies. While the basic steps of a pipeline are fairly universal -definition of brain regions, construction of an interaction matrix, and supervised learning-studies in the literature show many methodological variants (Table A2). Recommendations on methods that perform well can increase practitioner's productivity and limit vibration effects that risk undermining the reliability of biomarkers [START_REF] Varoquaux | Cross-validation failure: small sample sizes lead to large error bars[END_REF]. A challenge to such recommendations is the heterogeneity of prediction settings, for instance across different acquisition centers or clinical questions.

Here, we investigate methodological choices across 6 databases covering different clinical questions and behavioral task. We systematically compare commonly used functional connectome-based prediction methods. We find that some trends emerge, despite a large variance due to variability across subjects -visible across the folds-and across cohorts and clinical questions. Non-sparse linear models, such as logistic regression, appear as a good default choice of classifier. The lack of success of sparse approaches suggests that the discriminant signal is distributed across the functional connectome for the tasks we study. The tangent-space parametrization of functional connectomes brings improvements to prediction accuracy. With regards to nodes of the functional connectomes, defining them from rest-fMRI data gives slight benefits in prediction. Linear decomposition methods, such as dictionary learning or ICA, are good approaches to define these nodes from the rest-fMRI data at hand. Unlike clustering methods base on "hard" assignment, they provide a soft assignment to regions, enabling to capture a form of uncertainty in the definition of regions. Alternatively, the MODL15 (Mensch et al., 2016a) or BASC [START_REF] Bellec | Multi-level bootstrap analysis of stable clusters in resting-state fMRI[END_REF] atlases provide good readily-available nodes that simplify the process and alleviate computational cost. The good analytic performance of pre-computed atlases is promising and calls for further study. Establishing standard atlases brings significant computational benefits, as the definition of regions and the extraction of signal is the most computation-intensive part of the pipeline -in particular when performed inside a nested cross-validation loop. We found that using around 100 networks (corresponding to 150 regions) was sufficient for good prediction, though for many region-definition approaches a finer resolution did not hurt average prediction accuracy but only increased variance.

Overall, these results are consistent with the practice of the field. Preliminary comparisons in [START_REF] Abraham | Deriving reproducible biomarkers from multi-site resting-state data: An autism-based example[END_REF] on a single cohort revealed similar trends though ICA had performed poorly while here, with more systematic benchmarking, it appears to be a good solution. ICA has been used to define functional parcellations or nodes of functional connectomes by many groups [START_REF] Kiviniemi | Functional segmentation of the brain cortex using high model order group PICA[END_REF][START_REF] Rashid | Dynamic connectivity states estimated from resting fMRI identify differences among schizophrenia, bipolar disorder, and healthy control subjects[END_REF][START_REF] Smith | A positive-negative mode of population covariation links brain connectivity, demographics and behavior[END_REF][START_REF] Miller | Multimodal population brain imaging in the UK biobank prospective epidemiological study[END_REF]. More generally, it is well recognized that the nodes should be defined to match functional networks [START_REF] Smith | Network modelling methods for fMRI[END_REF]. Logistic regression, or the closely-related SVM, is the go-to classifier for many. Tangent-space parametrization of the connectivity matrix is more exotic, probably due to the mathematical complexity of its original presentation. However, it is gaining traction outside of methods studies [START_REF] Colclough | The heritability of multi-modal connectivity in human brain activity[END_REF][START_REF] Ng | Distinct alterations in parkinson's medicationstate and disease-state connectivity[END_REF] and is simple to implement, as summarized in Appendix A.

To enable comparison across different cohorts, we focused on 2-class classification problems. However, the results in terms of regions definition and connectivity parametrization should extend to other supervised learning settings, such as regression -e.g. for age prediction [START_REF] Liem | Predicting brain-age from multimodal imaging data captures cognitive impairment[END_REF]-multi-output approaches as with Canonical Correlation Analysis popular in large-scale population imaging settings [START_REF] Smith | A positive-negative mode of population covariation links brain connectivity, demographics and behavior[END_REF][START_REF] Miller | Multimodal population brain imaging in the UK biobank prospective epidemiological study[END_REF] for dimensional approaches to psychology.

Limitations and Challenges. The main limitation of our study is probably that we had to make choices and focus on the most popular methods. Indeed, to study systematically methods avoiding overfit requires computationalintense nested cross-validation (where the nesting is required to set the methods' internal parameters). In particular, we did not investigate Total-Variation constrained dictionary learning (TV-MSDL, [START_REF] Abraham | Extracting brain regions from rest fMRI with totalvariation constrained dictionary learning[END_REF]). This approach defines regions by imposing spatial structure in a linear-decomposition model. In a previous study, we found it promising [START_REF] Abraham | Deriving reproducible biomarkers from multi-site resting-state data: An autism-based example[END_REF], but it entailed too large of a computational cost for this multicohort study. Another important class of methods that this study did not investigate are biomarkers based on graph-theoretical approaches. Indeed, we benchmarked variants of a specific pipeline -region definition, followed by construction of a connectivity matrix, and supervised learning on it. Graph-theoretical approaches are an additional step to add to this pipeline. A full study of all options with this additional step would result in a combinatorial explosion of pipelines and prohibitive computational costs. We hope that the good choices of regions for edge-level models outlined in this study is also a good one for graph-theoretical approaches and that further studies can focus on exploring only a subset of the options covered here.

With evolving techniques, characteristics of data change, and optimal choices may evolve. However, the consistency of results on HCP suggest that our conclusions apply to high-quality datasets using state of the art techniques. A potential concern is the low accuracy for markers of drug abuse in subjects from ACPI datasets, possibly because the number of subjects is small or because ADHD status confounds drug-abuse predictions even after regressing it out. Nevertheless, our pipelines achieved similar accuracy as reported in a previous study on the same data [START_REF] Meszlényi | A model for classification based on the functional connectivity pattern dynamics of the brain[END_REF]. Finally, the analysis performed here can only outline trends across datasets. Indeed, the study does not establish that a pipeline choice strictly dominates others (see notes on statistical analysis

Step Recommendation

1: region extraction Functional regions, eg Dictionary learning or ICA 2: connectivity matrix Tangent-space embedding 3: supervised learning Non-sparse linear model, eg logistic regression or SVM in Appendix J), but it gives expected improvements. In term of expected improvement, the choice of classifier is the most important, as going from a poor to a good choice can improve the AUC by more than .1. Both choice of region and choice of parametrization bring smaller expected improvements.

Conclusion

Predictive models on rest-fMRI bring the promise of robust and reliable biomarkers: given new brain imaging data, they should give accurate predictions of clinics or behavior [START_REF] Woo | Building better biomarkers: brain models in translational neuroimaging[END_REF]. The framework of the functional connectomes grounds well the analysis of rest-fMRI; yet instantiating it still calls for many arbitrary choices.

Our study reveals trends that can provide good defaults to practitioners, summarized on Table 3: regions defined from functional data, for instance with ICA or dictionary learning as in the pre-computed MODL atlas, representing connectivity with the tangent embedding of covariance matrices, and using a non-sparse linear model, such as a logistic regression. In particular, good defaults can limit the combinatorial explosion of analytic pipelines, which decreases the computational cost of running a study and makes its conclusion more robust statistically. Yet, as it is well known in machine learning [START_REF] Wolpert | The lack of a priori distinctions between learning algorithms[END_REF], there cannot be a one-size-fits-all solution to data analysis: optimal choices will differ on datasets with very different properties from the datasets studied here. As can be seen, tangent-space parametrized connections are interpretable and positions in between correlation and partial correlation in terms of connectivity differences. We show the matrices estimated using timeseries extracted with pre-computed MODL dict. learning atlas of n=64.

Bristol-Myers

Appendix A. Computing the covariance tangentspace

Most of the methods that we study are readily-available in several computing environments, including Matlab and Python with a variety of well-maintained implementations. However, the only library that provides the tangent-space parametrization of covariance matrices is the Nilearn Python library16 . To facilitate reproducing our analysis in different environments, we describe here how to compute this parametrization with a few simple formulas. The computation is made of two step: First a group average covariance matrix Σ⋆ is computed from the covariances of the training subjects: {Σi, i ∈ T rain}. Second, it is used to transform covariance matrices, in the train set or the test set.

Computing the group average. As with any analysis based on covariance or correlation matrices, it is preferable to compute individual covariances from time series with an estimator that ensures well-conditioned matrices. The [START_REF] Ledoit | A well-conditioned estimator for largedimensional covariance matrices[END_REF] estimator is a good default choice [START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF][START_REF] Brier | Partial covariance based functional connectivity computation using ledoit-wolf covariance regularization[END_REF].

Strictly speaking, the group average should be computed according to the geometry of covariance matrices (Varoquaux et al., 2010a;[START_REF] Pennec | A Riemannian framework for tensor computing[END_REF]. This is a Frechet mean, which is computed by minimizing a cost function for instance using algorithm 3 of [START_REF] Fletcher | Riemannian geometry for the statistical analysis of diffusion tensor data[END_REF]. A simpler approach relies on using the Euclidean mean, which we found to give almost the same predictive performance. In this case, the formula of the mean is the standard one:

Euclidean mean ∶ Σ⋆ = 1 ntrain i∈Train Σi (A.1)
Transforming covariance matrices. Given the group reference covariance matrix Σ⋆, covariance matrices are transformed in the tangent-space representation by whitening them as follows (Varoquaux et al., 2010a). Computations are easily written with eigenvalues decompositions 17 : given a subject's covariance matrix Σi, 1. Compute the whitened matrix Σi

= Σ -1 2 ⋆ Σi Σ -1 2 ⋆ : Σi ← U T ∆ -1 2 U Σi U T ∆ -1 2 U (A.2)
where U T ∆ U = Σ⋆ by eigen-value decomposition, and operations on the diagonal matrix ∆ are element-wise operation applied to the diagonal. 2. Compute the matrix logarithm logm Σi:

logm( Σi) = ŨT log( ∆i ) Ũ (A.3)
where Σi = ŨT ∆i Ũ and the logarithm is applied to the diagonal elements of ∆i. Finally, the resulting matrix is turned to a vector and its entries are used as a features for the classifier.

The motivation from these transformations arises from the fact that covariance matrices -or correlations matrices-form a specific manifold of the R p×p matrices. Their structure is broken by standard additive arithmetic's: the difference of two covariances may create a matrix that does not correspond to the covariance matrix of a signal. Optimal statistical analysis calls for following the structure of the manifold [START_REF] Pennec | A Riemannian framework for tensor computing[END_REF]. The tangent-space parametrization is a simple way to approximate this structure by Euclidean geometry, in which standard additions and subtractions can be used (Varoquaux et al., 2010a).

With regards to statistical analysis, the structure of covariance matrices appears as constraints, or dependencies, between the coefficients of the matrix. As a result, these coefficients alone form a poor representation for second-level statistical analysis. The tangent-space approximation yields a parametrization of the problem in which features are i.i.d. (Varoquaux et al., 2010a). Such a parametrization is optimal for statistical learning. In addition, as discussed in Varoquaux et al. (2010a), this parametrization also gives good edge-level tests for instance see Figure A1. Hence, the weight vectors of the classifiers can be interpreted as edge-level weights. method, the distribution of the number of regions regularly increases with increasing dimensionality: the number of regions is roughly proportional to the number of networks.

For the optimal choice of dimensionality (dim = 80 for Group ICA and DictLearn and dim = 120 for K-Means as studied in Figure 6), the average number of regions lies on average around 150.

The number of regions extracted from ICA network sometimes displays an ill-controlled behavior, for instance for dim = 80, 100, 120 on COBRE, where a small number of folds lead to 300 or more regions. We believe that these high number of regions are extracted from an noisy ICA map with little structure as shown in Figure A8. Dictionary learning, which has a criteria on sparsity of of the maps, does not create such unstructured maps, and therefore does not suffer from the same problem.

In general, it can be expected that increasing the dimensionality of an ICA typically split large-scale networks into smaller regions. Given these noise components, the picture is slightly more complex. Indeed, as shown by Figure A7, high-dimensional ICA also typically extracts more noise components.

Table A1: 5th percentile, median and 95th percentile of accuracy scores in AUC over cross-validation folds (n = 100) for all six rest-fMRI datasets including HCP. Accuracy scores reported correspond to optimal choices in functional connectivity prediction pipeline as shown on Figure 7: precomputed atlas defined using massive online dictionary learning (MODL), connectivity matrices parametrized by their tangent-space representation, and an 2 -regularized logistic regression as a classifier.

parameter that sets the amount of sparsity. A good regularization parameter typically depends on the amount of time points available and the number of nodes. For each dataset and choice of regions, we used an inner-loop optimization to set the regularization parameter: we test two parameters values, 0.5 and 0.2, where 0 corresponds to a non regularized covariance and 1 or above to fully sparse covariances. For all datasets, 0.5 was the best trade-off as using 0.2 gave ill-conditioned results. Note that parameter selection for the graphical lasso on a wide and varied dataset is challenging as the graphical lasso runs into convergence problems when covering a variety of regularization parameters on covariance matrices with different properties. This happens very seldom, but with 100 folds on 1 500 subjects, and a dozen different atlases, the problem makes automatic parameter selection difficult. The convergence problem is well understood theoretically: the algorithm is a primal-dual algorithm and errors can accumulate between the primal and the dual solution [START_REF] Mazumder | The graphical lasso: New insights and alternatives[END_REF].

Figure A11 shows the impact of graph lasso on connectivity parametrizations. Overall, we observed no improvements in the prediction results based on graph lasso with respect to Ledoit-Wolf.

Figure 1 :

 1 Figure 1: Functional connectomeprediction pipeline with three main steps: 1) definition of brain regions (ROIs) from rest-fMRI images or using already defined reference atlases, 2) quantifying functional interactions from time series signals extracted from these ROIs and 3) comparisons of functional interactions across subjects using supervised learning.

Figure 2 :

 2 Figure 2: Brain regions extracted with ICA, DictLearn, KMeans, and Ward For ICA and dictionary learning, the dimensionality is of 80 and 60 resting-state networks -which are then broken up into more regions -yielding 150 regions, and 120 for KMeans and Ward clustering. Colors are arbitrary.

  Figure3summarizes the performances of classifiers on prediction scores for all rest-fMRI datasets. The results display a certain amount of variance across folds and datasets (i.e., prediction targets). However, they show that non-sparse ( 2 -regularized) linear classifiers perform better, with a slight lead for logistic-2 . Using non-linear classifiers does not appear useful; neither does sparsity. The results in Figure3are conditional on a good choice for the other steps of the pipeline. The marginal performances of the different choices of classifiers -i.e. cconsidering all other choices in the pipeline-are shown in Figure A2. They show similar trends, leading to preferring 2 -regularized linear classifiers.

Figure 6 :

 6 Figure6: Impact of the number of regions in atlases on prediction accuracy. The figure shows the distribution of the relative accuracy AUC scores across methods on the five rest-fMRI datasets, as a function of the number of regions. Horizontal bars (black) represent the median of the relative scores for the given number of regions. The chosen dimensionality for each method is indicated by a red arrow and was selected as the one with lowest variance in the error, and a median above zero.

Figure A1 :

 A1 Figure A1: Difference between mean of MCI and AD group connectivity matrices: We show the connectivity matrices from the ADNI dataset computed on samples diagnosed as Mild Cognitive Impairment (MCI) and Alzheimer's disease (AD).As can be seen, tangent-space parametrized connections are interpretable and positions in between correlation and partial correlation in terms of connectivity differences. We show the matrices estimated using timeseries extracted with pre-computed MODL dict. learning atlas of n=64.

Figure A2 :Figure A3 :

 A2A3 Figure A2: Impact of classifiers on prediction accuracy:Marginal distribution of the relative prediction accuracy of all classification choices for all rest-fMRI datasets. The results are obtained covering all the choices for the remaining steps i.e., atlases and connectivity parametrizations. Non-sparse linear models perform well and 2 -regularized logistic regression appears as the best choice.

Figure A6 :

 A6 FigureA6: Impact of the number of regions in data-driven atlases on prediction accuracy. For each method, the distribution of relative prediction accuracy (AUC) is displayed as a function of the number of regions, across rest-fMRI datasets. The horizontal bar (black) represents the median of the relatively mean scores for each dimensionality. The whiskers of each data point represent the 95% confidence interval for a given dataset across the folds. The optimal choices (red arrow) are selected as the minimal variance score with a median close to the maximum. The prediction scores are obtained for the optimal pipeline setting, involving tangent-space parameterization and logistic regression-2 classifier.

Table 1 :

 1 Datasets and prediction tasks, as well as the number of subjects in each group. COBRE -142 subjects, ADNI -136 subjects, ADNIDOD -167 subjects, ACPI -126 subjects, ABIDE -866 subjects, HCP -443. IQ represents fluid intelligence; 788 subjects had an IQ score in the HCP900 release. The acquisition parameters of each dataset are summarized in TableA3.and excluded from the analysis if they have severe scanner artifacts or head movements with amplitude larger than 2mm. Since pre-processed rest-fMRI subjects from ABIDE and ACPI are available, we choose images pre-

	10 , Gaussian
	spatial smoothing (F W HM = 5mm). The SPM based
	preprocessing pipeline is implemented through pyprepro-
	cess 11 -Python scripts relying on Nipype interface (Gor-
	golewski et al., 2011). All subjects were visually inspected
	6 www.adni-info.org
	7 www.adni-info.org/DOD.html

8 http://fcon_1000.projects.nitrc.org/indi/ACPI/html/

9

We perform some additional experiments on the Human Connectome Project (HCP) data, to assess that our experimental results still hold when using high-quality datasets like HCP.

10 www.fil.ion.ucl.ac.uk/spm/ 11 https://github.com/neurospin/pypreprocess processed using C-PAC pipeline

(Craddock et al.)

, without global signal regression. For ACPI, we choose linearly registered images using (Advanced Normalization Tools) ANTS and without motion scrubbing and no global signal regression. For already available preprocessed rest-fMRI subjects, we select the protocols such that it matches with the standard protocol we use. We have not done any additional preprocessing steps on ABIDE and ACPI.

Table 2 :

 2 5th percentile, median and 95th percentile of accuracy scores in AUC over cross-validation folds (n = 100) for all five rest-fMRI datasets.

Table 3 :

 3 Recommendations for rest-fMRI based prediction pipeline.

  Squibb Company; CereSpir, Inc.; Cogstate; Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Re-The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer's Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.

	search and Development, LLC.; Johnson and Johnson
	Pharmaceutical Research and Development LLC.; Lumos-
	ity; Lundbeck; Merck and Co., Inc.; Meso Scale Diagnos-
	tics, LLC.; NeuroRx Research; Neurotrack Technologies;
	Novartis Pharmaceuticals Corporation; Pfizer Inc.; Pira-
	mal Imaging; Servier; Takeda Pharmaceutical Company;
	and Transition Therapeutics.

We used the term ball rather a sphere. From a mathematical standpoint, A "ball" is the inside of a sphere.

Pre-computed sparse dictionaries with the MODL approach ofMensch et al. (2016a) are available from https://team.inria.fr/ parietal/files/2018/10/MODL_rois.zip

We also investigated higher dimensionality(150, 200 and 300) on some of the datasets, but could not do a systematic study above 300 because of high computational costs. These preliminary results showed no improvements in prediction accuracy compared to lower dimensionalities. It should be noted that the resulting components typically encompass several brain regions, which explains the dimension difference.

We also included Lasso as another choice of classifier in the pipeline. We observed significantly low prediction performance.

cobre.mrn.org

Note that, these curves are shown for the optimal choices found above: an 2 -penalized logistic regression as a classifier, and tangentspace parametrization to clarify the interpretation.

To ensure a correct nested cross-validation and avoid circularity (overfitting), data-driven region-extraction methods must be run on each fold, hence several hundred time for each pipeline configuration.

https://team.inria.fr/parietal/files/2018/10/MODL_rois. zip

Chen, G.,Ward, B.D., Xie, C., Li, W., Wu, Z., Jones, J.L., Franczak, M., Antuono, P., Li, S.J., 2011. Classification of alzheimer disease, mild cognitive impairment, and normal cognitive status with large-scale network analysis based on resting-state functional MR imaging. Radiology 259, 213-221.

http://nilearn.github.io/

http://preprocessed-connectomes-project.org/abide/

http://fcon_1000.projects.nitrc.org/indi/ACPI/html/

Acknowledgments

This work is funded by NiConnect project (ANR-11-BINF-0004 NiConnect) and the CATI project. We also thank to the open-source data community and preprocessed data initiatives for giving access to rest-fMRI datasets. Data collection and sharing for this project was funded by the Alzheimer's Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: Ab-bVie, Alzheimer's Association; Alzheimer's Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.; Biogen;

Data and code availability

The original data can be retrieved from: • COBRE http://cobre.mrn.org/ • ADNI & ADNIDOD www.adni-info.org • ACPI http://fcon_1000.projects.nitrc.org/indi/ ACPI/html/ • ABIDE http://preprocessed-connectomes-project. org/abide/index.html • HCP https://www.humanconnectome. org/study/hcp-young-adult/document/ 900-subjects-data-release Code to reproduce the experiments can be retrieved from: https://github.com/KamalakerDadi/ benchmark_rsfMRI_prediction

To facilitate reproduction, we also share time series, though without labels to comply with original data usage conditions, on https://github.com/ KamalakerDadi/benchmark_rsfMRI_prediction/blob/ master/downloader.py

Ethics statement

No experiments on living beings were performed for this study. Hence IRB approval was not necessary. The data that we used were acquired in original studies that had received approval by the original institution's IRB. All data were used accordingly to respective usage guidelines.

Appendix B. Time-series signals extraction

In this appendix, we give more details on time-series extraction, to complement subsection 2.3. After defining brain ROIs, we extract a representative time-series for each ROI in each subject. For atlases composed of non-overlapping ROIs 17 All covariance matrices are symmetric definite positive, and wellconditioned if estimated with the Ledoit-Wolf approach as can be seen in Figure 2 (bottom row), we simply compute the weighted average of the fMRI time series signals over all voxels within that specific region. For fuzzy overlapping ROIs, such as the atlases driven by CanICA and DictLearn as shown in Figure 2 (top row), we use ordinary least squares regression to unmix the signal in each voxel as the additive decomposition of signals over several overlapping ROIs. This is the same procedure as in [START_REF] Abraham | Deriving reproducible biomarkers from multi-site resting-state data: An autism-based example[END_REF]. Let Y ∈ R n×p be the subject-specific signals, written as p voxels by n timepoints, and V ∈ R k×p the atlas of k maps supported on p voxels. We estimate U ∈ R n×k , the set of time series for each ROI, using:

At the signal-extraction level, we regress out confounds or non-neural information [START_REF] Varoquaux | Learning and comparing functional connectomes across subjects[END_REF]. As confounding time-series we use: 10 CompCor [START_REF] Behzadi | A component based noise correction method (compcor) for BOLD and perfusion based fMRI[END_REF] on the whole brain and 6 motion related. We remove motion-related signal only for COBRE, ADNI and ADNIDOD as they are provided as raw data. We have not done any additional preprocessing steps on already preprocessed public datasets like ABIDE 18 , ACPI 19 . The signal of each region is also then normalized, detrended and bandpass-filtered between 0.01 and 0.1Hz. All these steps are done with Nilearn v0.3.

Investigating filtering choices.

At signal extraction level, we perform additional experiments to assess the impact of filtering strategies (low-pass filtering, global signal regression) on prediction accuracy. Overall, we observe no significant differences between filtering strategies -low-pass filter and no global signal regression, low-pass filter and global signal regression, no low-pass filter and no global signal regression. See Appendix F for complete comparisons.

Appendix C. Comparing each step marginally on the others

The figures in the main part of the paper summarize the impact of one modeling choice in the pipeline conditionally on nearly-optimal choices for the other steps. Here we compare the modeling choices at each step marginally on all other choices, i.e. considering all results and including well performing and poorly performing pipelines. This approach studies each step of the method independently from the other steps. The results are overall similar to performing the conditional analysis, however the variance is larger, as the plots pool together pipelines that perform well and pipelines that perform poorly.

Appendix C.1. Step:3 Choice of classifiers

Figure A2 shows the relative impact of classifier choices on prediction accuracy. All 2 regularized classifiers are performing markedly better than other considered classifiers, with a logistic regression as the best performer.

Appendix C.2. Step 2: Choices of connectivity parametrizations

Figure A3 shows the relative impact of connectivity parametrization on prediction accuracy when considering all choices for the other pipeline steps.

The tangent-space parametrization performs better than correlation or partial correlation and gives less variance.

Appendix C.3. Step:1 Choices of region-definition methods

Figure A4 shows the relative impact of region-definition choices on prediction accuracy for all choices in remaining steps. Atlases which are functionally derived lead to good performances. Linear-decomposition methods, appear as the best Figure A4: Impact of regions-definition method on prediction accuracy: Marginal distribution of relative prediction accuracy per region-definition approach across all rest-fMRI datasets. This is obtained while considering all pipelining options uniformly in all other steps. Among all pre-defined atlases, BASC is best. Among datadriven based atlases, the best choice is Online dictionary learning.

Appendix D. Additional experiments on regiondefinition methods

Here we give additional results related to step 1 of the pipeline: defining nodes, formed of brain regions or brain networks. An important parameter to chose is the optimal dimensionality dim i.e., how many networks are needed to predict from the rest-fMRI images. Another choice is whether these networks should be broken up into simply connected regionswith a region-extraction step-or whether distributed networks can be readily used. These two parameters may be important in the comparison of brain-region definition, in particular for data-driven approaches such as linear decomposition or clustering methods. In our study, we found that atlases learned using linear decomposition methods give a good prediction.

Appendix D.1. Varying dimensionality without region extraction

While the results in the main part of the manuscript are presented after regions extraction, Figure A6 studies the optimal dimensionality without region extraction. The optimal choices are not far from what we have observed with region extraction (Figure 6). As shown with a red arrow, the dimensions shown good prediction impact are: ICA and DictLearn -80 networks, K-Means and Ward -120 clusters, 122 networks with BASC. As dimensionality goes higher, we observe that the variance in prediction accuracy increases.

Appendix D.2. Defining nodes with regions or networks

One the optimal dimensionality chosen with and without region extraction, we compare for each method whether region-definition with regions or networks gives best prediction. Figure A5 summarizes the results but with no obvious clear-cut conclusions. The figure shows the distribution of differences in prediction scores between regions-based approaches and network-based approaches. There is a very slight tendency Spurious maps: Group ICA decomposition to favor region-based approaches, but the trend is not significant.

Appendix E. Regions extracted for different network-definition methods

Data-driven method tend to naturally extract networks rather than regions: ICA and dictionary learning give distributed networks, while KMeans gives clusters made of different connected components. Only Ward clustering readily gives ROIs as it has connectivity constraints.

We use a region-extraction procedure to go from networks to regions (Abraham et al., 2014a). Figure A7 shows the distribution of the number of regions extracted from networks obtained with different approaches for increasing dimensionality. The computational cost is too high to learn spatial maps for each split using Group ICA and Online Dictionary Learning. Hence we show outcomes for dimensions up to 120. For each

Appendix F. Experiments on filtering time-series

The results in the main part of the paper (eg Figure 3, Figure 4, Figure 5) as well as the marginal distribution figures shown on Appendix C, are established with low-pass filter of the time series and without global signal regression.

Here, we investigate the impact of filtering strategies on prediction accuracy. We compare three filtering schemes: lowpass filtering and no global signal regression, low-pass filtering and global signal regression, no low-pass filtering and no global signal regression. Figure A9 shows the outcome for all of these filtering combinations: there is no significant differences on prediction accuracy across all rs-fMRI datasets. : Impact of covariance estimator, for the different connectivity parametrizations: Marginal distribution of relative prediction scores per connectivity method across covariance structures for all rest-fMRI datasets. This is obtained by considering the good choices in dimensionality as studied in the main figures. Overall, using empirical covariance instead of Ledoit-Wolf shrinkage yield minor difference in prediction perfomances; note that partial correlation and tangent parametrization require the use of shrinkage. Graph lasso tends to perform lower compared with other estimators.

Appendix G. Tangent-space parametrization: Euclidean versus Geometric mean

We also compare two variants of tangent-space parametrization:

using the geometric (Frechet) mean and using the simple Euclidean mean for Σ ⋆ (Appendix A). Figure A10 highlights the difference between the relative prediction scores for each dataset. We compare between the optimal choices in brain atlas methods i.e. GroupICA and DictLearn.

Geometric mean based parametrization gives a slightly better prediction accuracy than Euclidean mean based parametrization.

Euclidean mean gives a reduced computational cost, but the Geometric mean is better justified in theory.

Appendix H.1. Empirical covariance: unregularized

We considered empirical covariance as our choice of unregularized covariance estimator. Using empirical covariance in the pipelines did not improve the predictions. We found little difference when compared with Ledoit-Wolf using full correlation as can be seen from Figure A11.

The empirical covariances are ill-conditioned matrices and not invertible, hence they cannot be used for partial correlations or tangent parametrization. Ledoit-Wolf shrinkage or sparse inverse covariance estimators overcome such limitations.

Appendix H.2. GraphLasso: 1 -regularized estimator Sparse inverse covariance is effective at recovering brain connectivity [START_REF] Smith | Network modelling methods for fMRI[END_REF][START_REF] Varoquaux | Brain covariance selection: better individual functional connectivity models using population prior[END_REF]. We use the graphical lasso [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] to estimate sparse covariance matrices, in order to study full correlation, partial correlation and tangent space parametrization in comparison with 2-regularized covariance estimator (Ledoit-Wolf)

Appendix I. A review of current practices in functional connectome-based classification

Table A2 summarizes about the list of methods used for prediction studies in diverse psychiatric diseases.

Appendix J. A note of statistical analysis of crossvalidation

Cross-validation cannot easily be used to reject null hypotheses when comparing analytic choices as the multiple perfold values that it gives are not independent (as discussed in appendices of [START_REF] Varoquaux | Cross-validation failure: small sample sizes lead to large error bars[END_REF]). Rather, they are resampling estimates, and therefore give a posterior on the prediction accuracy: the prediction accuracy is a property of the model fitted on the data at hand -see sec 8.4 [START_REF] Hastie | The elements of statistical learning[END_REF] for the link between resampling and Bayesian statistics. Comparing distribution of cross-validation accuracy thus cannot establish frequentist p-values -that the differences observed between models are not due to chance-but it can give posterior predictive distribution, and therefore the expected improvement of one model compared to another on new data.


Annette Luce

Y Saibi

Isabelle Séverin

Julien Boudon

Nadine Millot

Marie-Christine Chagnon

In vitro hazard assessment of nanoparticles developed for biomedical applications

à 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.

Nanotechnology is a growing sector in industry for twenty years and nanomaterials are currently used in numerous industrials applications. Their interest is based on their dimensions and their shape that give them very particular technological and biological properties. However, the development of these nanoparticles, their industrial preparation and integration into various products already involve an initial increased of human exposure (via inhalation, percutaneous and oral route). Another exposure pathway exists and is represented by drug injections. The risks associated with these new technologies and new products are still unclear. After a characterization step (zeta potential, TEM, TGA, DLS, XPS, IR measurements …), different toxicological endpoints were analyzed in vitro. However, bioassay protocols need to be adapted to nanoparticles. Cytotoxicity was assessed with the RNA synthesis inhibition assay, which is an early very sensitive sublethal bioassay. Oxidative stress was detected using DCF-DA assay for a sensitive and rapid quantitation of oxygen-reactive species. For genotoxicity, the comet assay performed in presence of Formamido-Pyrimidine-Glycosylase (FPG) detects low levels of DNA damage relative to DNA oxidative damage. Finally a short-term cell transformation assay (CTA) was performed, using the Bhas 42 cell line to distinguish between tumor initiators and promotors properties. Nanoparticles (SPIONs and titanate nanotubes (TiONts)) didn't induce any cytotoxic, genotoxic and oxidative response after a 24h treatment with concentrations up to 100 µg/mL in the human HepG2 cell line. However, TiONts induced a positive dose dependent response in the CTA assay at the initiation step (1-10 µg/mL). ROS production induced by SPIONs and TiONts at high concentrations.

DCF-DA assay

 Short term assay to determine the carcinogenicity of nanoparticles (90% of carcinogenic substances detected 7 ).

Cell Transformation Assay

Principle: Observations of morphological change of Bhas 42 cells in foci. A significant increase in the number of foci shows a potential carcinogenicity (initiator and/or promotor effect) 8 .  TiONts induced initiation and then a potential carcinogenicity, mutagenic assay is running in the laboratory using mammalian cells.

Negative foci

 Concerning the equivocal promotion data observed with nanoparticles;

experiments are currently performed to check this response.

Transmission Electron Microscopy

Principle: Fluorometric assay using the cell permeant 1',7'-dichlorodihydrofluorescein diacetate (H 2 DCFDA). Fluorescence is proportional to ROS production.

Principle: Detect single and double DNA strand breaks and alkali-labile sites after an electrophoresis in alkaline condition (pH = 13) and DNA staining with a fluorescent intercalating agent (propidium iodide) 5 .

  ab , Y. Saibi ab , I. Séverin b , J. Boudon a , N. Millot a , M. C. Chagnon b

Figure 1 :Figure 2 :

 12 Figure 1 : Evaluation of genotoxicity in HepG2 cells incubated 24 hours in the presence of different concentrations of SPIONs. 3 independent experiments were performed and the results were averaged and expressed as score 6 .



  Simple assay to detect ROS production by HepG2 cell metabolism.

  Figure 4 : Fluorescence induced by ROS production after a 24h exposure to SPIONs -NH 2.

Figure 8 :

 8 Figure 7 : Number of foci /well induced after a 3 day exposure with TiONts in the growth phase in the initiation assay.

Figure 6 :Figure 9 :

 69 Figure 6 : Initiation and promotion protocols in the Bhas 42 cell transformation assay.

Figure 5 :9

 5 Figure 5 : TEM images of HepG2 cells incubated 4 hours in the presence of titanate nanotubes (a-b) and iron oxide nanoparticles (c-d) to a concentration of 100 µg/mL. Figures above show that the TiONts and the SPIONs have a vesicular localization.