Gaël Varoquaux

Reddy Raamana

Denis A Engemann

Andrés A Hoyos-Idrobo

Yannick A Schwartz

Bertrand A Thirion

Pradeep Reddy

Assessing and tuning brain decoders: cross-validation, caveats, and guidelines

Keywords: cross-validation, decoding, fMRI, model selection, sparse, bagging, MVPA

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: decoding needs model evaluation

Decoding, ie predicting behavior or phenotypes from brain images or signals, has become a central tool in neuroimage data processing [START_REF] Haynes | A primer on pattern-based approaches to fMRI: Principles, pitfalls, and perspectives[END_REF][START_REF] Haynes | Decoding mental states from brain activity in humans[END_REF][START_REF] Kamitani | Decoding the visual and subjective contents of the human brain[END_REF][START_REF] Norman | Beyond mind-reading: multi-voxel pattern analysis of fMRI data[END_REF][START_REF] Varoquaux | How machine learning is shaping cognitive neuroimaging[END_REF][START_REF] Yarkoni | Choosing prediction over explanation in psychology: Lessons from machine learning[END_REF]. In clinical applications, prediction opens the door to diagnosis or prognosis [START_REF] Demirci | A review of challenges in the use of fMRI for disease classification/characterization and a projection pursuit application from a multi-site fMRI schizophrenia study[END_REF][START_REF] Fu | Pattern classification of sad facial processing: toward the development of neurobiological markers in depression[END_REF][START_REF] Mouro-Miranda | Classifying brain states and determining the discriminating activation patterns: Support vector machine on functional MRI data[END_REF]. To study cognition, successful prediction is seen as evidence of a link between observed behavior and a brain region [START_REF] Haxby | Distributed and overlapping representations of faces and objects in ventral temporal cortex[END_REF] or a small fraction of the image [START_REF] Kriegeskorte | Information-based functional brain mapping[END_REF]. Decoding power can test if an encoding model describes well multiple facets of stimuli [START_REF] Mitchell | Predicting human brain activity associated with the meanings of nouns[END_REF][START_REF] Naselaris | Encoding and decoding in fMRI[END_REF]. Prediction can be used to establish what specific brain functions are implied by observed activations [START_REF] Poldrack | Decoding the large-scale structure of brain function by classifying mental states across individuals[END_REF][START_REF] Schwartz | Mapping cognitive ontologies to and from the brain[END_REF]. All these applications rely on measuring the predictive power of a decoder.

Assessing predictive power is difficult as it calls for characterizing the decoder on prospective data, rather than on the data at hand. Another challenge is that the decoder must often choose between many different estimates that give rise to the same prediction error on the data, when there are more features (voxels) than samples (brain images, trials, or subjects). For this choice, it relies on some form of regularization, that embodies a prior on the solution [START_REF] Hastie | The elements of statistical learning[END_REF]. The amount of regularization is a parameter of the decoder that may require tuning. Choosing a decoder, or setting appropriately its internal parameters, condition or the health status of subjects. The accuracy, or predictive power, of this model is the expected error on the prediction, formally: accuracy = E E(y pred , y ground truth )

where E is a measure of the error, most often1 the fraction of instances for which y pred ≠ y ground truth . Importantly, in equation ( 1), E denotes the expectation, ie the average error that the model would make on infinite amount of data generated from the same experimental process.

In decoding settings, the investigator has access to labeled data, ie brain images for which the variable to predict, y, is known. These data are used to train the model, fitting the model parameters, and to estimate its predictive power. However, the same observations cannot be used for both. Indeed, it is much easier to find the correct labels for brain images that have been seen by the decoder than for unknown images2 . The challenge is to measure the ability to generalize to new data.

The standard approach to measure predictive power is cross-validation: the available data is split into a train set, used to train the model, and a test set, unseen by the model during training and used to compute a prediction error (figure 1). Chapter 7 of [START_REF] Hastie | The elements of statistical learning[END_REF] contains a reference on statistical aspects of cross-validation. Below, we detail important considerations in neuroimaging.

Independence of train and test sets. Cross-validation relies on independence between the train and test sets. With time-series, as in fMRI, the autocorrelation of brain signals and the temporal structure of the confounds imply that a time separation is needed to give truly independent observations. In addition, to give a meaningful estimate of prediction power, the test set should contain new samples displaying all confounding uncontrolled sources of variability. For instance, in multi-session data, it is harder to predict on a new session than to leave out part of each session and use these samples as a test set. However, generalization to new sessions is useful to capture actual invariant information. Similarly, for multi-subject data, predictions on new subjects give results that hold at the population level. However, a confound such as movement may correlate with the diagnostic status predicted. In such a case the amount of movement should be balanced between train and test set. Sufficient test data. Large test sets are necessary to obtain sufficient power for the prediction error for each split of cross-validation. As the amount of data is limited, there is a balance to strike between achieving such large test sets the data is split multiple times into a train set, used to train the model, and a test set, used to compute predictive power.

Test set Train set

Full data and keeping enough training data to reach a good fit with the decoder. Indeed, theoretical results show that crossvalidation has a negative bias on small dataset [2, sec. [START_REF] Breiman | Submodel selection and evaluation in regression. the x-random case[END_REF].1] as it involves fitting models on a fraction of the data. On the other hand, large test sets decrease the variance of the estimated accuracy [2, sec.5.2]. A good cross-validation strategy balances these two opposite effects. Neuroimaging papers often use leave one out cross-validation, leaving out a single sample at each split. While this provides ample data for training, it maximizes test-set variance and does not yield stable estimates of predictive accuracy 3 . From a decision-theory standpoint, it is preferable to leave out 10% to 20% of the data, as in 10-fold cross-validation [18, chap. 7.12] [5, 27]. Finally, it is also beneficial to increase the number of splits while keeping a given ratio between train and test set size. For this purpose k-fold can be replaced by strategies relying on repeated random splits of the data (sometimes called repeated learning-testing4 [START_REF] Arlot | A survey of cross-validation procedures for model selection[END_REF] or ShuffleSplit [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]). As discussed above, such splits should be consistent with the dependence structure across the observations (using eg a LabelShuffleSplit), or the training set could be stratified to avoid class imbalance [START_REF] Raamana | Thickness network features for prognostic applications in dementia[END_REF]. In neuroimaging, good strategies often involve leaving out sessions or subjects.

Hyper-parameter selection

A necessary evil: one size does not fit all. In standard statistics, fitting a simple model on abundant data can be done without the tricky choice of a meta-parameter: all model parameters are estimated from the data, for instance with a maximum-likelihood criterion. However, in high-dimensional settings, when the number of model parameters is much larger than the sample size, some form of regularization is needed. Indeed, adjusting model parameters to best fit the data without restriction leads to overfit, ie fitting noise [18, chap. 7]. Some form of regularization or prior is then necessary to restrict model complexity, e.g. with low-dimensional PCA in discriminant analysis [START_REF] Chen | Exploring predictive and reproducible modeling with the single-subject FIAC dataset[END_REF], or by selecting a small number of voxels with a sparse penalty [START_REF] Carroll | Prediction and interpretation of distributed neural activity with sparse models[END_REF][START_REF] Yamashita | Sparse estimation automatically selects voxels relevant for the decoding of fMRI activity patterns[END_REF]. If too much regularization is imposed, the ensuing models are too constrained by the prior, they underfit, ie they do not exploit the full richness of the data. Both underfitting and overfitting are detrimental to predictive power and to the estimation of model weights, the decoder maps. Choosing the amount of regularization is a typical bias-variance problem: erring on the side of variance leads to overfit, while too much bias leads to underfit. In general, the best tradeoff is a data-specific choice, governed by the statistical power of the prediction task: the amount of data and the signal-to-noise ratio.

Nested cross-validation. Choosing the right amount of regularization can improve the predictive power of a decoder and controls the appearance of the weight maps. The most common approach to set it is to use crossvalidation to measure predictive power for various choices of regularization and to retain the value that maximizes predictive power. Importantly, with such a procedure, the amount of regularization becomes a parameter adjusted on data, and thus the predictive performance measured in the corresponding cross-validation loop is not a reliable assessment of the predictive performance of the model. The standard procedure is then to refit the model on the available data, and test its predictive performance on new data, called a validation set. Given a finite amount of data, a nested cross-validation procedure can be employed: the data are repeatedly split in a validation set and a decoding set to perform decoding. The decoding set itself is split in multiple train and test sets with the same validation set, forming an inner "nested" cross-validation loop used to set the regularization hyper-parameter, while the external loop varying the validation set is used to measure prediction performance -see figure 2.

Model averaging. Choosing the best model in a family of good models is hard. One option is to average the predictions of a set of suitable models [44, chap. 35], [START_REF] Churchill | Comparing within-subject classification and regularization methods in fMRI for large and small sample sizes[END_REF][START_REF] Hoyos-Idrobo | Improving sparse recovery on structured images with bagged clustering[END_REF][START_REF] Kuncheva | Classifier ensembles for fMRI data analysis: an experiment[END_REF] -see [18, chap. 8] for a description outside of neuroimaging. A simple version of this idea is bagging [START_REF] Breiman | Bagging predictors[END_REF]: using bootstrap, random resamplings of the data, to generate many train sets and corresponding models, the predictions of which are then averaged. The benefit of these approaches is that if the errors of each model are sufficiently independent, they average out: the average model performs better and displays much less variance. With linear models often used as decoders in neuroimaging, model averaging is appealing as it boils down to averaging weight maps.

To benefit from the stabilizing effect of model averaging in parameter tuning, we can use a variant of both crossvalidation and model averaging 5 . In a standard crossvalidation procedure, we repeatedly split the data in train and test set and for each split, compute the test error for a grid of hyper-parameter values. However, instead of selecting the hyper-parameter value that minimizes the mean test error across the different splits, we select for each split the model that minimizes the corresponding test error and average these models across splits.

Model selection for neuroimaging decoders

Decoding in neuroimaging faces specific modelselection challenges. The main challenge is probably the scarcity of data relative to their dimensionality, typically hundreds of observations 6 . Another important aspect of decoding is that, beyond predictive power, interpreting model weights is relevant.

Common decoders and their regularization. Both to prefer simpler models and to facilitate interpretation, linear models are ubiquitous in decoding. In fact, their weights form the common brain maps for visual interpretation.

The classifier used most often in fMRI is the support vector machine (SVM) [START_REF] Chen | Exploring predictive and reproducible modeling with the single-subject FIAC dataset[END_REF][START_REF] Laconte | Support vector machines for temporal classification of block design fMRI data[END_REF][START_REF] Mouro-Miranda | Classifying brain states and determining the discriminating activation patterns: Support vector machine on functional MRI data[END_REF]. However, logistic regressions (Log-Reg) are also often used [START_REF] Chen | Exploring predictive and reproducible modeling with the single-subject FIAC dataset[END_REF][START_REF] Rasmussen | Model sparsity and brain pattern interpretation of classification models in neuroimaging[END_REF][START_REF] Ryali | Sparse logistic regression for whole-brain classification of fMRI data[END_REF][START_REF] Varoquaux | Small-sample brain mapping: sparse recovery on spatially correlated designs with randomization and clustering[END_REF][START_REF] Yamashita | Sparse estimation automatically selects voxels relevant for the decoding of fMRI activity patterns[END_REF]. Both of these classifiers learn a linear model by minimizing the sum of a loss L -a data-fit term-and a penalty p -the regularizing energy term that favors simpler models:

ŵ = argmin w L(w) + 1 C p(w) L = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ SVM logistic
where C is the regularization parameter that controls the bias-variance tradeoff: small C means strong regularization. The SVM and logistic regression model differ only by the loss used. For the SVM the loss is a hinge loss: flat and exactly zero for well-classified samples and with a misclassification cost increasing linearly with distance to the literature. It is commonly used in other areas of machine learning, for instance to set parameters in bagged models such as trees, by monitoring the out-of-bag error (eg in the scikit-learn library [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]). 6 While in imaging neuroscience, hundreds of observations seems acceptably large, it is markedly below common sample sizes in machine learning. Indeed, data analysis in brain imaging has historically been driven by very simple models while machine learning has tackled rich models since its inception.

Large C

Small C C=1.

Log-reg 1 decision boundary. For the logistic regression, it is a logistic loss, which is a soft, exponentially-decreasing, version of the hinge [START_REF] Hastie | The elements of statistical learning[END_REF]. By far the most common regularization is the 2 penalty. Indeed, the common form of SVM uses 2 regularization, which we will denote SVM-2 . Combined with the large zero region of the hinge loss, strong 2 penalty implies that SVMs build their decision functions by combining a small number of training images (see figure 3). Logistic regression is similar: the loss has no flat region, and thus every sample is used, but some very weakly.

C=10 2 C=10 5 C=10 -5 SVM 2 C=1 L R C=10 5
Another frequent form of penalty, 1 , imposes sparsity on the weights: a strong regularization means that the weight maps w are mostly comprised of zero voxels (see Fig. 4).

Parameter-tuning in neuroimaging. In neuroimaging, many publications do not discuss their choice of decoder hyper-parameters; while others state that they use the default value, eg C = 1 for SVMs. Standard machine learning practice advocates setting the parameters by nested crossvalidation [START_REF] Hastie | The elements of statistical learning[END_REF]. For non sparse, 2 -penalized models, the amount of regularization often does not have a strong influence on the weight maps of the decoder (see figure 4). Indeed, regularization in these models changes the fraction of input maps supporting the hyperplane (see 3). As activation maps for the same condition often have similar aspects, this fraction impacts weakly decoders' maps.

For sparse models, using the 1 penalty, sparsity is often seen as a means to select relevant voxels for prediction [START_REF] Carroll | Prediction and interpretation of distributed neural activity with sparse models[END_REF][START_REF] Ryali | Sparse logistic regression for whole-brain classification of fMRI data[END_REF]. In this case, the amount of regularization has a very visible consequence on weight maps and voxel selection (see figure 4). Neuroimaging studies often set it by cross-validation [START_REF] Carroll | Prediction and interpretation of distributed neural activity with sparse models[END_REF], though very seldom nested (exceptions comprise [START_REF] Churchill | Comparing within-subject classification and regularization methods in fMRI for large and small sample sizes[END_REF][START_REF] Varoquaux | Small-sample brain mapping: sparse recovery on spatially correlated designs with randomization and clustering[END_REF]). Voxel selection by 1 penalty on brain maps is unstable because neighboring voxels respond similarly and 1 estimators will choose somewhat randomly few of these correlated features [START_REF] Rondina | Stabilitybased multivariate mapping using scors[END_REF][START_REF] Varoquaux | Small-sample brain mapping: sparse recovery on spatially correlated designs with randomization and clustering[END_REF]. Hence various strategies combining sparse models are used in neuroimaging to improve decoding performance and stability. Averaging weight maps across cross-validation folds [START_REF] Hoyos-Idrobo | Improving sparse recovery on structured images with bagged clustering[END_REF][START_REF] Varoquaux | Small-sample brain mapping: sparse recovery on spatially correlated designs with randomization and clustering[END_REF], as described above, is interesting, as it stays in the realm of linear models. Relatedly, [START_REF] Grosenick | Interpretable whole-brain prediction analysis with graphnet[END_REF] report the median of weight maps, thought it does not correspond to weights in a predictive model. Consensus between sparse models over data perturbations gives theoretically better feature selection [START_REF] Meinshausen | Stability selection[END_REF]. In fMRI, it has been used to screen voxels before fitting linear models [START_REF] Rondina | Stabilitybased multivariate mapping using scors[END_REF][START_REF] Varoquaux | Small-sample brain mapping: sparse recovery on spatially correlated designs with randomization and clustering[END_REF] or to interpret selected voxels [START_REF] Yamashita | Sparse estimation automatically selects voxels relevant for the decoding of fMRI activity patterns[END_REF].

For model selection in neuroimaging, prediction performance is not the only relevant metric and some control over the estimated model weights is also important. For this purpose, [START_REF] Laconte | The evaluation of preprocessing choices in single-subject bold fMRI using npairs performance metrics[END_REF][START_REF] Rasmussen | Model sparsity and brain pattern interpretation of classification models in neuroimaging[END_REF][START_REF] Strother | The quantitative evaluation of functional neuroimaging experiments: the npairs data analysis framework[END_REF] advocate using a tradeoff between prediction performance and stability of decoder maps. Stability is a proxy for estimation error on these maps, a quantity that is not accessible without knowing the ground truth. While very useful it gives only indirect information on estimation error: it does not control whether all the predictive brain regions were found, nor whether all regions found are predictive. Indeed, a decoder choosing its maps independently from the data would be very stable, though likely with poor prediction performance. Hence the challenge is in finding a good prediction-stability tradeoff [START_REF] Rasmussen | Model sparsity and brain pattern interpretation of classification models in neuroimaging[END_REF][START_REF] Strother | Stability and reproducibility in fMRI analysis[END_REF].

Empirical studies: cross-validation at work

Here we highlight practical aspects of cross-validation in brain decoding with simple experiments. We first demonstrate the variability of prediction estimates on MRI, MEG, and simulated data. We then explore how to tune decoders parameters.

Experiments on real neuroimaging data

A variety of decoding datasets. To achieve reliable empirical conclusions, it is important to consider a large number of different neuroimaging studies. We investigate crossvalidation in a large number of 2-class classification problems, from 7 different fMRI datasets (an exhaustive list can be found in Appendix E). We decode visual stimuli within subject (across sessions) in the classic Haxby dataset [START_REF] Haxby | Distributed and overlapping representations of faces and objects in ventral temporal cortex[END_REF], and across subjects using data from [START_REF] Duncan | Consistency and variability in functional localisers[END_REF]. We discriminate across subjects i) affective content with data from [START_REF] Wager | Neural mechanisms of emotion regulation: evidence for two independent prefrontal-subcortical pathways[END_REF], ii) visual from narrative with data from [START_REF] Moran | Social-cognitive deficits in normal aging[END_REF], iii) famous, familiar, and scrambled faces from a visual-presentations dataset [START_REF] Henson | Face repetition effects in implicit and explicit memory tests as measured by fMRI[END_REF], and iv) left and right saccades in data from [START_REF] Knops | Recruitment of an area involved in eye movements during mental arithmetic[END_REF]. We also use a non-published dataset, ds009 from openfMRI [START_REF] Poldrack | Toward open sharing of task-based fMRI data: the OpenfMRI project[END_REF]. All the across-subject predictions are performed on trial-by-trial response (Z-score maps) computed in a first-level GLM. Finally, beyond fMRI, we perform prediction of gender from VBM maps using the OASIS data [START_REF] Marcus | Open access series of imaging studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF]. Note that all these tasks cover many different settings, range from easy discriminations to hard ones, and (regarding fMRI) recruit very different systems with different effect size and variability. The number of observations available to the decoder varies between 80 (40 per class) and 400, with balanced classes.

The results and figures reported below are for all these datasets. We use more inter-subject than intra-subject datasets. However 15 classification tasks out of 31 are intra-subject (see Tab. A1). In addition, when decoding is performed intra-subject, each subject gives rise to a crossvalidation. Thus in our cross-validation study, 82% of the data points are for intra-subject settings.

All MR data but [START_REF] Knops | Recruitment of an area involved in eye movements during mental arithmetic[END_REF] are openly available from openfMRI [START_REF] Poldrack | Toward open sharing of task-based fMRI data: the OpenfMRI project[END_REF] or OASIS [START_REF] Marcus | Open access series of imaging studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF]. Standard preprocessing and first-level analysis were applied using SPM, Nipype and Nipy (details in Appendix F.1). All MR data were variance-normalized 7 and spatially-smoothed at 6 mm FWHM for fMRI data and 2 mm FWHM for VBM data.

MEG data. Beyond MR data, we assess cross-validation strategies for decoding of event-related dynamics in neurophysiological data. We analyze magneteoencephalography (MEG) data from a working-memory experiment made available by the Human Connectome Project [START_REF] Larson-Prior | Adding dynamics to the human connectome project with MEG[END_REF]. We perform intra-subject decoding in 52 subjects with two runs, using a temporal window on the sensor signals (as in [START_REF] Sitt | Large scale screening of neural signatures of consciousness in patients in a vegetative or minimally conscious state[END_REF]). Here, each run serves as validation set for the other run. We consider two-class decoding problems, focusing on either the image content (faces vs tools) or the functional role in the working memory task (target vs low-level and high-level distractors). This yields in total four classification analyzes per subject. For each trial, the feature set is a time window constrained to 50 ms before and 300 ms after event onset, emphasizing visual components. We use the cleaned single-trial outputs from the HCP "tmegpreproc" pipeline. MEG data analysis was performed with the MNE-Python software [START_REF] Gramfort | MEG and EEG data analysis with MNE-Python[END_REF][START_REF] Gramfort | MNE software for processing MEG and EEG data[END_REF]. Full details on the analysis are given in Appendix F.3.

Experimental setup. Our experiments make use of nested cross-validation for an accurate measure of prediction power. As in figure 2, we repeatedly split the data in a validation set and a decoding set passed on to the decoding procedure (including parameter-tuning for experiments in 3.3 and 3.4). To get a good measure of predictive power, we choose large validation sets of 50% of the data, respecting the sample dependence structure (leaving out subjects, 7 Division of each time series voxel/MEG sensor by its standard deviation or sessions). We use 10 different validation sets that each contribute a data point in results.

We follow standard decoding practice in fMRI [START_REF] Pereira | Machine learning classifiers and fMRI: a tutorial overview[END_REF]. We use univariate feature selection on the training set to select the strongest 20% of voxels and train a decoder on the selected features. As a choice of decoder, we explore classic linear models: SVM and logistic regression with 1 and 2 penalty 8 . We use scikit-learn for all decoders [START_REF] Abraham | Machine learning for neuroimaging with scikit-learn[END_REF][START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF].

In a first experiment, we compare decoder performance estimated by cross-validation on the decoding set, with performance measured on the validation set. In a second experiment, we investigate the use of cross-validation to tune the model's regularization parameter, either using the standard refitting approach, or averaging as described in section 2.2, as well as using the default C = 1 choice of parameter, and a value of C = 1000.

Results: cross-validation to assess predictive power

Reliability of the cross-validation measure. Considering that prediction error on the large left-out validation set is a good estimate of predictive power, we use it to assess the quality of the estimate given by the nested crossvalidation loop. Figure 5 shows the prediction error measured by cross-validation as a function of the validationset error across all datasets and validation splits. It reveals a small negative bias: as predicted by the theory, 8 Similar decoders adding a regularization that captures spatiotemporal correlations among the voxels are well suited for neuroimaging [START_REF] Gramfort | Identifying predictive regions from fMRI with TV-L1 prior[END_REF][START_REF] Grosenick | Interpretable whole-brain prediction analysis with graphnet[END_REF][START_REF] Karahanoglu | Total activation: fMRI deconvolution through spatiotemporal regularization[END_REF][START_REF] Michel | Total variation regularization for fMRI-based prediction of behavior[END_REF]. Also, random forests, based on model averaging discussed above, have been used in fMRI [START_REF] Kuncheva | Classifier ensembles for fMRI data analysis: an experiment[END_REF][START_REF] Langs | Detecting stable distributed patterns of brain activation using gini contrast[END_REF]. However, this review focuses on the most common practice. Indeed, these decoders entail computational costs that are intractable given the number of models fit in the experiments.

Cross-validation

Difference cross-validation is pessimistic compared to a model fit on the complete decoding set. However, models that perform poorly are often reported with a better performance by cross-validation. Additionally, cross-validation estimates display a large variance: there is a scatter between estimates in the nested cross-validation loop and in the validation set.

Different cross-validation strategies. Figure 6 summarizes the discrepancy between prediction accuracy measured by cross validation and on the validation set for different cross-validation strategies: leaving one sample out, leaving one block of data out -where blocks are the natural units of the experiment, sessions or subjects-and random splits leaving out 20% of the blocks of data with 3, 10, and 50 repetitions.

Ideally, a good cross-validation strategy would minimize this discrepancy. We find that leave-one-sample out is very optimistic in within-subject settings. This is expected, as samples are highly correlated. When leaving out blocks of data that minimize dependency between train and test set, the bias mostly disappears. The remaining discrepancy appears mostly as variance in the estimates of prediction accuracy. For repeated random splits of the data, the larger the number of splits, the smaller the variance. Performing 10 to 50 splits with 20% of the data blocks left out gives a better estimation than leaving successively each blocks out, at a fraction of the computing cost if the number of blocks is large. While intra and inter subject settings do not differ strongly when leaving out blocks of data, intra-subject settings display a larger variance of estimation as well as a slight negative bias. These are likely due to non-stationarity in the timeseries, eg scanner drift or loss of vigilance. In inter-subject settings, heterogeneities may hinder prediction [START_REF] Strother | Stability and reproducibility in fMRI analysis[END_REF], yet a cross-validation strategy with multiple subjects in the test set will yield a good estimate of prediction accuracy 9 .

Other modalities: MEG and simulations. We run the experiments on the MEG decoding tasks and the simulations. We generate simple simulated data that mimic brain imaging to better understand trends and limitations of cross-validation. Briefly, we generate data with 2 classes in 100 dimensions with Gaussian noise temporally autocorrelated and varying the separation between the class centers (more details in Appendix A). We run the experiments on a decoding set of 200 samples.

The results, displayed in figure 7, reproduce the trends observed on MR data. As the simulated data is temporally auto-correlated, leave-one-sample-out is strongly optimistic. Detailed analysis varying the separability of the classes (Appendix A) shows that cross-validation tends to be pessimistic for high-accuracy situations, but optimistic when prediction is low. For MEG decoding, the leave-oneout procedure is on trials, and thus does not suffer from correlations between samples. Cross-validation is slightly pessimistic and display a large variance, most likely because of inhomogeneities across samples. In both situations, leaving blocks of data out with many splits (e.g. 50) gives best results.

Results on cross-validation for parameter tuning

We now evaluate cross-validation as a way of setting decoder hyperparameters.

Tuning curves: opening the black box. Figure 8 is a didactic view on the parameter-selection problem: it gives, for varying values of the meta-parameter C, the crossvalidated error and the validation error for a given split of validation data 10 . The validation error is computed on a large sample size on left out data, hence it is a good estimate of the generalization error of the decoder. Note that the parameter-tuning procedure does not have access to this information. The discrepancy between the tuning curve, computed with cross-validation on the data available to the decoder, and the validation curve, is an indication of the uncertainty on the cross-validated estimate of prediction power. Test-set error curves of individual splits of the nested cross-validation loop show plateaus and a discrete behavior. Indeed, each individual test set contains dozens of observations. The small combinatorials limit the accuracy of error estimates.

Figure 8 also shows that non-sparse -2 -penalizedmodels are not very sensitive to the choice of the regularization parameter C: the tuning curves display a wide plateau 11 . However, for sparse models ( 1 models), the maximum of the tuning curve is a more narrow peak, particularly so for SVM. A narrow peak in a tuning curve implies that a choice of optimal parameter may not alway carry over to the validation set.

Impact of parameter tuning on prediction accuracy. Crossvalidation is often used to select regularization hyperparameters, eg to control the amount of sparsity. On figure 9, we compare the various strategies: refitting with the best parameters selected by nested cross-validation, averaging the best models in the nested cross-validation, or simply using either the default value of C or a large one, given that tuning curves can plateau for large C.

For non-sparse models, the figure shows that tuning the hyper-parameter by nested cross validation does not lead in general to better prediction performance than a default choice of hyper-parameter. Detailed investigations (figure A3) show that these conclusions hold well across all tasks, though refitting after nested cross-validation is beneficial for good prediction accuracies, ie when there is either a large signal-to-noise ratio or many samples.

For sparse models, the picture is slightly different. Indeed, high values of C lead to poor performancepresumably as the models are overly sparse-, while using default value C = 1, refitting or averaging models tuned by cross-validation all perform well. Investigating how these compromises vary as a function of model accuracy (figure A3) reveals that for difficult decoding situations (low prediction) it is preferable to use the default C = 1, while in good decoding situations, in is beneficial to tune C by nested cross-validation and rely on model averaging, which tends to perform well and displays less variance.

Results: stability of model weights

Impact of parameter tuning on stability. The choice of regularization parameter also affects the stability of the weight maps of the classifier. Strongly regularized maps underfit, thus depending less on the train data, which may lead to increased stability. We measure stability of the decoder maps by computing their correlation across different choices of validation split for a given task.

Figure 10 summarizes the results on stability. For all models, sparse and non-sparse, model averaging does give more stable maps, followed by refitting after choosing parameters by nested cross-validation. Sparse models are much less stable than non-sparse ones [START_REF] Varoquaux | Small-sample brain mapping: sparse recovery on spatially correlated designs with randomization and clustering[END_REF].

Prediction power -stability tradeoff. The choice of decoder with the best predictive performance might not give the most stable weight maps, as seen by comparing figures 9 and 10. Figure 11 shows the prediction-stability tradeoff for different decoders and different parameter-tuning strategies. Overall, SVM and logistic-regression perform similarly and the dominant effect is that of the regularization: non-sparse, 2 -penalized, models are much more stable than sparse, 1 -penalized, models. For non-sparse models, averaging stands out as giving a large gain in stability albeit with a decrease of a couple of percent in prediction accuracy compared to using C = 1 or C = 1000, which gives good prediction and stability (figures 9 and 11). With sparse models, averaging offers a slight edge for stability and, for SVM performs also well in prediction. C = 1000 achieves low stability (figure 10), low prediction power (figure 9), and a poor tradeoff.

Appendix D.2 shows trends on datasets where the prediction is easy or not. For non-sparse models, averaging brings a larger gain in stability when prediction accuracy is large. Conversely, for sparse models, it is more beneficial to average in case of poor prediction accuracy.

Note that these experimental results are for common neuroimaging settings, with variance-normalization and univariate feature screening. [START_REF] Haxby | Distributed and overlapping representations of faces and objects in ventral temporal cortex[END_REF]. The thin colored lines are test scores for each of the internal cross-validation folds, the thick black line is the average of these test scores on all folds, and the thick dashed line is the score on left-out validation data. The vertical dashed line is the parameter selected on the inner cross-validation score.

Discussion and conclusion: lessons learned

Decoding seeks to establish a predictive link between brain images and behavioral or phenotypical variables. Prediction is intrinsically a notion related to new data, and therefore it is hard to measure. Cross-validation is the tool of choice to assess performance of a decoder and

C V + a v e

Sparse models

Figure 10: Stability of the weights: impact of the parametertuning strategy: for each strategy, difference to the mean stability of the model weights, where the stability is the correlation of the weights across validation splits. As the stability is a correlation, the unit is a different between correlation values. The reference is the mean stability across all models for a given prediction task.

to tune its parameters. The strength of cross-validation is that it relies on few assumptions and probes directly the ability to predict, unlike other model-selection procedures -eg based on information theoretic or Bayesian criteria. However, it is limited by the small sample sizes typically available in neuroimaging 12 .

An imprecise assessment of prediction. The imprecision on the estimation of decoder performance by crossvalidation is often underestimated. Empirical confidence intervals of cross-validated accuracy measures typically extend more than 10 points up and down (figure 5 and6). Experiments on MRI (anatomical and functional), MEG, and simulations consistently exhibit these large error bars due to data scarcity. Such limitations should be kept in mind for many MVPA practices that use predictive power as a form of hypothesis testing -eg searchlight [START_REF] Kriegeskorte | Information-based functional brain mapping[END_REF] or testing for generalization [START_REF] Knops | Recruitment of an area involved in eye movements during mental arithmetic[END_REF]-and it is recommended to use permutation to define the null hypothesis [START_REF] Kriegeskorte | Information-based functional brain mapping[END_REF]. In addition, in the light of cross-validation variance, methods publications should use several datasets to validate a new model.

Guidelines on cross-validation. Leave-one-out crossvalidation should be avoided, as it yields more variable results. Leaving out blocks of correlated observations, rather than individual observations, is crucial for nonbiased estimates. Relying on repeated random splits with 20% of the data enables better estimates with less computation by increasing the number of cross-validations without shrinking the size of the test set.

Parameter tuning. Selecting optimal parameters can improve prediction and change drastically the aspects of weight maps (Fig. 4).

However, our empirical study shows that for variancenormalized neuroimaging data, non-sparse decoders ( 2penalized) are only weakly sensitive to the choice of their parameter, particularly for the SVM. As a result, relying on the default value of the parameter often outperforms parameter tuning by nested cross-validation. Yet, such parameter tuning tends to improve the stability of the maps. For sparse decoders ( 1 -penalized), default parameters also give good prediction performance. However, parameter tuning with model averaging increases stability and can lead to better prediction. Note that it is often useful to variance normalize the data (see Appendix C).

Concluding remarks. Evaluating a decoder is hard. Crossvalidation should not be considered as a silver bullet. Neither should prediction performance be the only metric. To assess decoding accuracy, best practice is to use repeated learning-testing with 20% of the data left out, while keeping in mind the large variance of the procedure. Any parameter tuning should be performed in nested crossvalidation, to limit optimistic biases. Given the variance that arises from small samples, the choice of decoders and their parameters should be guided by several datasets.

Our extensive empirical validation (31 decoding tasks, with 8 datasets and almost 1 000 validation splits with nested cross-validation) shows that sparse models, in particular 1 SVM with model averaging, give better prediction but worst weight-maps stability than non-sparse classifiers. If stability of weight maps is important, non-sparse SVM with C = 1 appears to be a good choice. Further work calls for empirical studies of decoder performance with more datasets, to reveal factors of the dataset that could guide better the choice of a decoder for a given task. Results. Figure A2 summarizes the cross-validation measures for different values of separability. Beyond the effect of the cross-validation strategy observed on other figures, the effect of the separability, ie the true prediction accuracy is also visible. Setting aside the leave-one-sample-out cross-validation strategy, which is strongly biased by the correlations across the samples, we see that all strategies tend to be biased positively for low accuracy and negatively for high accuracy. This observation is in accordance with trends observed on figure 5.

Cross-validation

Appendix B. Comparing parameter-tuning strategies

Figure A3 shows pairwise comparisons of parameter-tuning strategies, in sparse and non-sparse situations, for the bestperforming options. In particular, it investigates when different strategies should be preferred. The trends are small. Yet, it appears that for low predictive power, setting C=1 in nonsparse models is preferable to cross-validation while for high predictive power, cross-validation is as efficient. This is consistent with results in figure 5 showing that cross-validation is more reliable to measure prediction error in situations with a good accuracy than in situations with a poor accuracy. Similar trends can be found when comparing to C=1000. For sparse models, model averaging can be preferable to refitting. We however find that for low prediction accuracy it is favorable to use C=1, in particular for logistic regression.

Sparse models

Figure A3: Comparing parameter-tuning strategies on prediction accuracy. Relating pairs of tuning strategies. Each dot corresponds to a given dataset, task, and validation split. The line is an indication of the tendency, using a lowess non-parametric local regression. The top row summarizes results for non-sparse models, SVM 2 and logistic regression 2 ; while the bottom row gives results for sparse models, SVM 1 and logistic regression 1 .

Note these figures show points related to different studies and classification tasks. The trends observed are fairly homogeneous and there are not regions of the diagram that stand out. Hence, the various conclusions on the comparison of decoding strategies are driven by all studies.

Appendix C. Results without variancenormalization

Results without variance normalization of the voxels are given in figure A4 for the correspondence between error measured in the inner cross-validation loop, figure A5 for the effect of the choice of a parameter-tuning strategy on the prediction performance, and figure A6 for the effect on the stability of the weights.

Cross-validation on non variance-normalized neuroimaging data is not more reliable than on variance-normalized data (figure A4). However, parameter tuning by nested cross-validation is more important than on variance-normalized data for good prediction (figure A5). This difference can be explained by the fact that variance normalizing makes dataset more comparable to each others, and thus a default value of parameters is more likely to work well.

In conclusion, variance normalizing the data can be important, in particular with non-sparse SVM. Table A1: The different datasets and tasks. We report the prediction performance on the validation test for parameter tuning by 10 random splits followed by refitting using the best parameter.

All analyses were based on the cleaned single-trial outputs obtained from the HCP "tmegpreproc" pipeline which provides cleaned segmented sensor space data. The MEG data that were recorded with a wholehead MAGNES 3600 (4D Neuroimaging, San Diego, CA) magnetometer system in a magnetically shielded room. Contamination by environmental magnetic fields was accounted for by computing the residual MEG signal from concomitant recordings of reference gradiometers and magnetometers located remotely from the main sensor array. Data were bandpass filtered between 1.3 and 150Hz using zero-phase forward and everse Butterworth filters. Notch filters were then applied at (59-61/119-121 Hz) to attenuate line noise artefacts. Data segments contaminated by remaining environmental or system artifacts were detected using a semi-automatic HCP pipeline that takes into account the local and global variation as well as the correlation structure of the data. Independent component analysis based on the FastICA algorithm was then used to estimate and suppress spatial patterns of cardiac and ocular artifacts. Artifact related components were identified in a semi-automatic fashion assisted by comparisons with concomitantly recorded electrocardiogram (ECG) and electrooculogram (EOG). These components were then projected out from the data. Depending on the classification of bad channels performed by the HCP pipelines, the data contained fewer than 248 sensors. For details on the HCP pipelines see Larson-Prior et al. [START_REF] Larson-Prior | Adding dynamics to the human connectome project with MEG[END_REF] and the HCP reference manual. The MEG data were accessed through the MNE-Python software [START_REF] Gramfort | MEG and EEG data analysis with MNE-Python[END_REF][START_REF] Gramfort | MNE software for processing MEG and EEG data[END_REF] and the MNE-HCP library .

Appendix G. Performance on each classification task

The prediction accuracy results presented in the various figures are differential effects removing the contribution of the dataset. In figure A9, we present for each decoding strategy the prediction accuracy on all datasets.

We can see that the variations of prediction accuracy from one decoding strategy to another are mostly reported across datasets: the various lines are roughly parallel. One notable a v e ra g in g re fi tt in g C = 1 C = 1 0 0 0 a v e ra g in g re fi tt in g C = 1 C = 1 0 0 0 a v e ra g in g re fi tt in g C = 1 C = 1 0 0 0 a v e ra g in g re fi tt in g C = 1 C = 1 0 0 0 exception is SVM 1 with C = 1000 for which some datasets show a strong decrease. Another, weaker, variation is the fact that 1 models tend to perform better on the Haxby dataset (our source of intra-subject classification tasks). This good performance of sparse models could be due to the intra-subject settings: sparse maps are less robust to inter-subject variability. However, the core messages of the paper relative to which parameter-tuning strategy to use are applicable to intra and inter-subject settings. For non-sparse models, using a large value of C without parameter tuning is an overall safe choice, and for sparse models, model averaging, refitting, and a choice of C = 1 do not offer a clear win, although model averaging is comparatively less variable.

Figure 1 :

 1 Figure 1: Cross-validation: the data is split multiple times into a train set, used to train the model, and a test set, used to compute predictive power.

Figure 2 :

 2 Figure 2: Nested cross-validation: two cross-validation loops are run one inside the other.

Figure 3 :

 3 Figure 3: Regularization with SVM-2 : blue and brown points are training samples of each class. The SVM learns a separating line between the two classes. In a weakly regularized setting (large C, this line is supported by few observations -called support vectors-, circled in black on the figure, while in a strongly-regularized case (small C), it is supported by the whole data.

C=10 - 5

 5

Figure 4 :

 4 Figure 4: Varying amount of regularization on the face vs house discrimination in the Haxby 2001 data[START_REF] Haxby | Distributed and overlapping representations of faces and objects in ventral temporal cortex[END_REF]. Left: with a log-reg 1 , more regularization (small C) induces sparsity. Right: with an SVM 2 , small C means that weight maps are a combination of a larger number of original images, although this has only a small visual impact on the corresponding brain maps.

Figure 5 :

 5 Figure 5: Prediction error: cross-validated versus validation set.Each point is a measure of predictor error in the inner crossvalidation loop (10 splits, leaving out 20%), or in the left-out validation set. The dark line is an indication of the tendency, using a lowess local regression. The two light lines indicate the boundaries above and below which fall 5% of the points. They are estimated using a Gaussian kernel density estimator, with a bandwidth set by the Scott method, and computing the CDF along the y direction.

Figure 7 :

 7 Figure 7: Cross-validation error: non-MRI modalities. Difference between accuracy measured by cross-validation and on the validation set, for MEG data and simulated data, with different crossvalidation strategies. Detailed simulation results in Appendix A.

1 Figure 8 :

 18 Figure 8: Tuning curves for SVM 2 , logistic regression 2 , SVM 1 , and logistic regression 1 , on the scissors / scramble discrimination for the Haxby dataset[START_REF] Haxby | Distributed and overlapping representations of faces and objects in ventral temporal cortex[END_REF]. The thin colored lines are test scores for each of the internal cross-validation folds, the thick black line is the average of these test scores on all folds, and the thick dashed line is the score on left-out validation data. The vertical dashed line is the parameter selected on the inner cross-validation score.

  r a g in g C V + r e fi tt in g C

Figure 9 :

 9 Figure 9: Prediction accuracy: impact of the parametertuning strategy. For each strategy, difference to the mean prediction accuracy in a given validation split.

Figure 11 :

 11 Figure 11: Prediction -stability tradeoff The figure summarize figures 9 and 10, reporting the stability of the weights, relative to the split's average, as a function of the delta in prediction accuracy. It provides an overall summary, with errors bar giving the first and last quartiles (more detailed figures in Appendix D)

Figure A4 :Figure A5 :Figure A6 :

 A4A5A6 Figure A4: Prediction error: cross-validated versus validation set. Each point is a measure of predictor error measure in the inner cross-validation loop, or in the left-out validation dataset for a model refit using the best parameters. The line is an indication of the tendency, using a lowess non-parametric local regression.

Figure A7 :

 A7 Figure A7: Prediction -stability tradeoff This figures gives the data points behind 11, reporting the stability of the weights, relative to the split's average, as a function of the delta in prediction accuracy. Each point is corresponds to a specific prediction task in our study.

Figure A8 :

 A8 Figure A8: Prediction -stability tradeoff The figure reports for each dataset and task the stability of the weights as a function of the prediction accuracy measured on the validation set, with the different choices of decoders and parameter-tuning strategy. The line is an indication of the tendency, using a lowess local regression. The stability of the weights is the correlation across validation splits.

Figure A9 :

 A9 FigureA9: Performance of each decoding strategy and each dataset. the plot is a "parallel coordinate plot": each lines denotes a dataset and the different x positions denote different decoding strategies.

  Cross-validation error: different strategies. Difference between accuracy measured by cross-validation and on the validation set, in intra and inter-subject settings, for different crossvalidation strategies: leave one sample out, leave one block of samples out (where the block is the natural unit of the experiment: subject or session), or random splits leaving out 20% of the blocks as test data, with 3, 10, or 50 random splits. For inter-subject settings, leave one sample out corresponds to leaving a session out. The box gives the quartiles, while the whiskers give the 5 and 95 percentiles.

				in accuracy measured
	strategy	by cross-validation and on validation set
	sample out Leave one		22%	+19% +3%	+43%
	subject/session Leave one		10% 21%	+10% +17%
	3 splits 20% leftout,		11% 24%	+11% +16%
	20% leftout, 10 splits 20% leftout, 50 splits		24% 23% 9% 9%	+14% +9% +13% +7%	Intra subject Inter subject
		40%	20% 10% 0% +10% +20%	+40%
	cross-validation < validation set	cross-validation > validation set
	Figure 6:

  FigureA2: Cross-validation measures on simulations. Prediction accuracy, as measured by cross-validation (box plots) and on a very large test set (vertical lines) for different separability on the simulated data and for different cross-validation strategies: leave one sample out, leave one block of samples out (where the block is the natural unit of the experiment: subject or session), or random splits leaving out 20% of the blocks as test data, with 3, 10, or 50 random splits. The box gives the quartiles, while the whiskers give the 5 and 95 percentiles. Note that here leave-one-block-out is similar of 10 splits of 10% of the data.

		Accuracy measured by cross-validation
	strategy Leave one sample out	0.05	for different separability 86% 94% Separability 0.1 0.2	97% 99%
	57%		72%			98% 100%
	Leave one block out		67%			87%
	20% leftout, 3 splits	55%	66%	76%		83%	89%	99%
	57%		73%	81%		100%
	20% leftout, 10 splits		68%				87%
		58%	71%		84%		99%
	20% leftout, 50 splits		69%			86%
	assymptotic prediction accuracy	62.9% sep=0.05	77.2% sep=0.1	84%	94.0% sep=0.2	98%

For multi-class problems, where there is more than

categories in y, or for unbalanced classes, a more elaborate choice is advisable, to distinguish misses and false detections for each class.2 A simple strategy that makes no errors on seen images is simply to store all these images during the training and, when asked to predict on an image, to look up the corresponding label in the store.

One simple aspect of the shortcomings of small test sets is that they produce unbalanced dataset, in particular leave-one-out for which there is only one class represented in the test set.

Also related is bootstrap CV, which may however duplicate samples inside the training set of the test set.

The combination of cross-validation and model averaging is not new (see eg[START_REF] Hoyos-Idrobo | Improving sparse recovery on structured images with bagged clustering[END_REF]), but it is seldom discussed in the neuroimaging

The probability of correct classification for each subject is also an interesting quantity, though it is not the same thing as the prediction accuracy measured by cross-validation [18, sec 7.12]. It can be computed by non-parametric approaches such as bootstrapping the train set[START_REF] Platt | Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods[END_REF], or using a posterior probability, as given by certain classifiers.

For the figure, we compute cross-validated error with a leaveone-session-out on the first 6 sessions of the scissor / scramble Haxby data, and use the last 6 sessions as a validation set.

[START_REF] Fu | Pattern classification of sad facial processing: toward the development of neurobiological markers in depression[END_REF] This plateau is due to the flat, or nearly flat, regions of their loss that renders them mostly dependent only on whether samples are well classified or not.

There is a trend in inter-subject analysis to acquire databases with a larger number of subjects, eg ADNI, HCP. Conclusions of our empirical study might not readily transfer to these settings.

https://github.com/nilearn/nilearn/graphs/contributors

https://github.com/scikit-learn/scikit-learn/graphs/ contributors

https://github.com/mne-tools/mne-python/graphs/ contributors

the values we explore for µ were chosen empirically to vary classification accuracy from 60% to 90%.

Wellcome Department of Cognitive Neurology, http://www.fil. ion.ucl.ac.uk/spm

Acknowledgments

This work was supported by the EU FP7/2007-2013 under grant agreement no. 604102 (HBP). Computing resource were provided by the NiConnect project (ANR-11-BINF-0004 NiConnect) and an Amazon Webservices Research Grant. The authors would like to thank the developers of nilearn 13 , scikit-learn 14 and MNE-Python 15 for continuous efforts in producing high-quality tools crucial for this work.

In addition, we acknowledge useful feedback from Russ Poldrack on the manuscript.

We generate data with samples from two classes, each described by a Gaussian of identity covariance in 100 dimensions. The classes are centered respectively on vectors (µ, . . . , µ) and (-µ, . . . , -µ) where µ is a parameter adjusted to control the separability of the classes. With larger µ the expected predictive accuracy would be higher. In addition, to mimic the time dependence in neuroimaging data we apply a Gaussian smoothing filter in the sample direction on the noise (σ = 2). Code to reproduce the simulations can be found on https: //github.com/GaelVaroquaux/cross_val_experiments.

We produce different datasets with predefined separability by varying 16 µ in (.05, .1, .2). Figure A1 shows two of these configurations.

Appendix A.2. Experiments: error varying separability

Unlike with a brain imaging datasets, simulations open the door to measuring the actual prediction performance of a classifier, and therefore comparing it to the cross-validation measure.

For this purpose, we generate a pseudo-experimental data with 200 train samples, and a separate very large test set, with 10 000 samples. The train samples correspond to the data available during a neuroimaging experiment, and we perform crossvalidation on these. We apply the decoder on the test set. The large number of test samples provides a good measure of prediction power of the decoder [START_REF] Arlot | A survey of cross-validation procedures for model selection[END_REF]. We use the same decoders as for brain-imaging data and repeat the whole procedure 100 times. For cross-validation strategies that rely on sample blocks -as with sessions-, we divide the data in 10 continuous blocks.

10% 5%

0% +5% 11 summarizes the effects of the decoding strategy on the prediction -stability tradeoff. On figure A7, we give the data points that underly this summary.

Appendix D.2. Prediction and stability interactions

Figure 11 captures the effects of the decoding strategy across all datasets. However, some classification tasks are easier or more stable than others.

We give an additional figure, figure A8, showing this interaction between classification performance and the best decoding strategy in terms of weight stability. The main factor of variation of the prediction accuracy is the choice of dataset, ie the difficulty of the prediction task.

Here again, we see that the most important choice is that of the penalty: logistic regression and SVM have overall the same behavior. In terms of stability of the weights, higher prediction accuracy does correspond to more stability, except for overly-penalized sparse model (C=1000). For non-sparse models, model averaging after cross-validation is particularly beneficial in good prediction situations.

Appendix E. Details on datasets used

Table A1 lists all the studies used in our experiments, as well as the specific prediction tasks. In the Haxby dataset [START_REF] Haxby | Distributed and overlapping representations of faces and objects in ventral temporal cortex[END_REF] we use various pairs of visual stimuli, with differing difficulty. We excluded pairs for which decoding was unsuccessful, such as scissors versus bottle).

Appendix F. Details on preprocessing

Appendix F.1. fMRI data Intra-subject prediction. For intra-subject prediction, we use the Haxby dataset [START_REF] Haxby | Distributed and overlapping representations of faces and objects in ventral temporal cortex[END_REF] as provided from the PyMVPA [START_REF] Hanke | PyMVPA: A python toolbox for multivariate pattern analysis of fmri data[END_REF] website -http://dev.pymvpa.org/datadb/haxby2001.html. Details of the preprocessing are not given in the original paper, beyond the fact that no spatial smoothing was performed. We have not performed additional preprocessing on top of this publicly-available dataset, aside from spatial smoothing with an isotropic Gaussian kernel, FWHM of 6 mm (nilearn 0.2, Python 2.7).

Inter-subject prediction. For inter-subject prediction, we use different datasets available on openfMRI [START_REF] Poldrack | Toward open sharing of task-based fMRI data: the OpenfMRI project[END_REF]. For all the datasets, we performed standard preprocessing with SPM8 17 : in the following order, slice-time correction, motion correction (realign), corregistration of mean EPI on subject's T1 image, and normalization to template space with unified segmentation on the T1 image. The preprocessing pipeline was orchestrated through the Nipype processing infrastructure [START_REF] Gorgolewski | Nipype: a flexible, lightweight and extensible neuroimaging data processing framework in python[END_REF]. For each subject, we then performed session-level GLM with a design according to the individual studies, as described in the openfMRI files, using Nipy (version 0.3, Python version 2.7) [START_REF] Millman | Analysis of functional magnetic resonance imaging in python[END_REF].

Appendix F.2. Structural MR data

For prediction from structural MR data, we perform Voxel Based Morphometry (VBM) on the Oasis dataset [START_REF] Marcus | Open access series of imaging studies (OASIS): cross-sectional MRI data in young, middle aged, nondemented, and demented older adults[END_REF]. We use SPM8 with the following steps: segmentation of the white matter/grey matter/CSF compartiments and estimation of the deformation fields with DARTEL [START_REF] Ashburner | Diffeomorphic registration using geodesic shooting and gauss-newton optimisation[END_REF]. The inputs for predictive models is the modulated grey-matter intensity. The corresponding maps can be downloaded with the dataset-downloading facilities of the nilearn software (function nilearn.datasets.fetch oasis vbm).

Appendix F.3. MEG data

The magneteoencephalography (MEG) data is from an Nback working-memory experiment made available by the Human Connectome Project [START_REF] Larson-Prior | Adding dynamics to the human connectome project with MEG[END_REF]. Data from 52 subjects and two runs was analyzed using a temporal window approach in which all magnetic fields sampled by the sensor array in a fixed time interval yield one variable set (see for example [START_REF] Sitt | Large scale screening of neural signatures of consciousness in patients in a vegetative or minimally conscious state[END_REF] for event related potentials in electroencephalography). Here, each of the two runs served as validation set for the other run. For consistency, two-class decoding problems were considered, focussing on either the image content (faces VS tools) or the functional role in the working memory task (target VS low-level and highlevel distractors). This yielded in total four classification analyses per subject. For each trial, the time window was then constrained to 50 millisecond before and 300 millisecond after event onset, emphasizing visual components.


domain: sde chim

text:
M Nihemaiti
email: maolida.nihemaiti@kaust.edu.sa

J Le Roux

J-P Croué

Formation of emerging disinfection byproducts by chlorination/chloramination of seawater impacted by algal organic matter

Keywords: Chlorine, chloramines, disinfection byproducts, algal organic matter, haloacetamides

The aim of this work was to study the formation of haloacetamides (HAcAms) and other DBPs during chlorination and chloramination of algal organic matter (AlOM). The HAcAms formation potentials of different precursors (amino acids, simulated algal blooms grown in the Red Sea) were evaluated. Experiments with simulated algal blooms were conducted in the presence of bromide ion (synthetic seawater containing 800 μg/L Br -) to assess the formation of brominated analogues of HAcAms in conditions close to the disinfection of real seawater. Chlorination produced more HAcAms than chloramination from real algae (Synecococcus sp.), thus indicating that the nitrogen of HAcAms comes predominantly from DON through the decarboxylation of amino acids rather than from NH2Cl. Dibrominated species of DBPs (i.e., DBAcAm, DBAA, DBAN) were the dominant species formed by both chlorination and chloramination of algal blooms samples. Chloramination of the amino acid asparagine produced an important amount of DCAcAm as compared to chlorination, indicating the existence of a specific reaction pathway.

INTRODUCTION

Prechlorination is used in seawater desalination processes to control biofouling in both thermal plants and membrane plants. Important amounts of disinfectant can be required to maintain a residual during algal blooms events. Many disinfection byproducts (DBPs) are formed during seawater pretreatment processes using chlorination, including the regulated trihalomethanes (THMs) and haloacetic acids (HAAs). Nitrogenous DBPs (N-DBPs) generally form in much smaller amounts than regulated DBPs, but have been a growing concern because of their greater health risk [START_REF] Plewa | Halonitromethane drinking water disinfection byproducts: Chemical characterization and mammalian cell cytotoxicity and genotoxicity[END_REF][START_REF] Muellner | Haloacetonitriles vs. regulated haloacetic acids: Are nitrogen-containing DBPs more toxic?[END_REF]. In vitro mamallian cell tests demonstrated that haloacetonitriles (HANs), halonitromethanes (HNMs) and haloacetamides (HAcAms) are more cytotoxic and genotoxic (up to 2 orders of magnitude) than non-nitrogenous THMs and HAAs [START_REF] Plewa | Comparative Mammalian Cell Toxicity of N-DBPs and C-DBPs. In Disinfection By-Products in Drinking Water[END_REF].

As an alternative disinfectant, monochloramine (NH2Cl) can be used in cooling systems to get a stable disinfectant residual. However, the use of NH2Cl can be a source of nitrogen and leads to the formation of N-DBPs. Moreover, the disinfection of seawater can produce important amounts of N-DBPs, especially during algal blooms, enriched in nitrogen-containing compounds (e.g., amino acids). Arid regions produce water by desalination of seawater or brackish water, characterized by elevated bromide and iodide ions contents. The presence of high concentrations of bromide (60 mg Br -/L in the Red Sea) and iodide (0.05 mg I -/L) ion in seawater favors the formation of brominated and iodinated byproducts that are often more toxic than their chlorinated analogues (Richardson et al., 2007, Richardson et al., 2008).

Haloacetamides are N-DBPs that exert high levels of toxicity compared to other DBPs, but their formation mechanisms are unclear. They are known to be products of the hydrolysis of HANs [START_REF] Reckhow | Formation and degradation of dichloroacetonitrile in drinking waters[END_REF], however they were reported to be formed independently by chlorination and chloramination of natural organic matter (NOM) [START_REF] Huang | Dichloroacetonitrile and dichloroacetamide can form independently during chlorination and chloramination of drinking waters, model organic matters, and wastewater effluents[END_REF]. Two reaction pathways have been proposed for the formation of haloacetamides during water disinfection. The first one is the decarboxylation pathway, producing haloacetonitriles that are then hydrolyzed into haloacetamides. In this case, the nitrogen atom of haloacetamides is coming from the precursor compound (i.e., amino acids, organic nitrogen) [START_REF] Huang | Dichloroacetonitrile and dichloroacetamide can form independently during chlorination and chloramination of drinking waters, model organic matters, and wastewater effluents[END_REF]. The second reaction pathway is the aldehyde pathway, where NH2Cl reacts with aldehydes to produce HANs and HAcAms [START_REF] Kimura | Chloroacetonitrile and N,2dichloroacetamide formation from the reaction of chloroacetaldehyde and monochloramine in water[END_REF]. In this case, the nitrogen of HAcAms is incorporated from NH2Cl. Literature about the formation mechanism of HAcAms and their precursors in NOM, wastewater, algal organic matter (AlOM) or seawater organic matter is limited. The HAcAm formation potential of some amino acids [START_REF] Chu | Precursors of dichloroacetamide, an emerging nitrogenous DBP formed during chlorination or chloramination[END_REF][START_REF] Huang | Dichloroacetonitrile and dichloroacetamide can form independently during chlorination and chloramination of drinking waters, model organic matters, and wastewater effluents[END_REF], humic acids, wastewater effluent and algal EPS [START_REF] Huang | Dichloroacetonitrile and dichloroacetamide can form independently during chlorination and chloramination of drinking waters, model organic matters, and wastewater effluents[END_REF] and culture bacteria [START_REF] Huang | Formation of haloacetonitriles and haloacetamides during chlorination of pure culture bacteria[END_REF] were reported but no studies were focused on seawater algal organic matter.

This study investigates the formation mechanisms of haloacetamides during chlorination and chloramination of algal organic matter (AlOM). The aim was to evaluate the importance of the decarboxylation and the aldehyde pathway from different precursors (amino acids, simulated algal blooms). Experiments with simulated algal blooms were conducted in the presence of bromide (synthetic seawater containing 60-65 mg/L Br -) to assess the formation of brominated analogues of haloacetamides.

MATERIALS AND METHODS

Materials

All experiments were conducted using deionized water (Milli-Q, Millipore) buffered with a mixture of sodium phosphate monobasic and sodium phosphate dibasic. Solution pH values were adjusted as needed using sodium hydroxide or sulfuric acid (0.1 N, Fisher Scientific). Fisher Scientific MTBE and ethyl acetate (>99%) were used without further purification. Amino acids (L-asparagine >98%, L-tyrosine >99%, L-tryptophan >99% and phenylalanine >99%) were used without further purification and were supplied through Sigma-Aldrich (Figure 1). Sodium hypochlorite (NaOCl, 5.65-6%, Fisher Scientific) and ammonium chloride (Acros Organics, 99.6%) were used as chlorination reagents. Ammonium-15 N-chloride was purchased from Sigma-Aldrich. Sodium thiosulfate (Fisher Scientific) was used to quench residual oxidant. A mixed standard containing haloacetonitriles (HANs), trichloronitromethane (TCNM) and haloketones (HKs) (EPA 551B Halogenated Volatiles Mix), a mixed standard containing 9 HAAs (EPA 552.2 Methyl Ester Calibration Mix) and surrogate standard decafluorobiphenyl (99%) were supplied from Supelco (Sigma-Aldrich). Chloro-, bromo-dichloro-and trichloroacetamide were obtained from Sigma-Aldrich. Other haloacetamides (HAcAms) were purchased from Cansyn Chem. Corp.

Table 1. Molecular structures of investigated amino acids L-asparagine L-tryptophan L-tyrosine Phenylalanine

Preparation and analysis of chlorine and chloramine Monochloramine (NH2Cl) stock solutions (50 mM) were prepared fresh daily by dissolving ammonium chloride (NH4Cl) in deionized water adjusted to pH = 8.5 with sodium hydroxide. Sodium hypochlorite (NaOCl) was then added slowly to the rapidly stirred solution, at a Cl:N molar ratio of at least 1:1.2 to avoid breakpoint chlorination resulting from local excess of hypochlorite. Adjusting the pH at 8.5 minimizes the disproportionation of NH2Cl to dichloramine (NHCl2), since NHCl2 forms at pH < 8 according to equilibrium 1:

2NH2Cl + H + = NHCl2 + NH4 + (1)
NH2Cl and NHCl2 were quantified in stock solutions by monitoring absorbances at their respective max (NH2Cl = 245 nm; NHCL2 = 295 nm). Residual oxidant in samples after experiments was analysed iodometrically (Standards Methods for the Examination of Water and Wastewater, 1995).

Experimental methods

All glassware used during these experiments was washed and baked at 500 °C for at least 5 h prior to use. Reactions were conducted in sealed 65 mL amber glass bottles at room temperature (20 °C). Chloramination experiments were conducted following the approach of [START_REF] Krasner | Development of DBP and Nitrosamine Formation Potential Tests for Treated Wastewater, Reclaimed Water, and Drinking Water[END_REF], using an excess of oxidant (15 mg/L as Cl2) and a reaction time of 72 hours at pH 8.

Simulated algal blooms grown directly in the Red Sea (mesocosms) were employed as algal source (main species Synecococcus sp.). Aqueous solutions were prepared by dilution of seawater samples (simulated algal blooms) with synthetic seawater (Grassoff et al., 1976) to the desired DOC concentration (1.5 mg C/L). Amino acids solutions were prepared by dissolving a pre-determined amount of amino acids in Milli-Q water containing 10 mM carbonate buffer. 275 µL of preformed monochloramine was then added. Each series of experiments included a blank (buffered Milli-Q water) and each sample was triplicated. After 72h of contact time, 5 mL of samples were used for residual chlorine analysis and 3 x 20 mL were used for DBPs analysis. Molar yields were calculated based on the molar concentrations of amino acids (mol DBP formed × 100 / mol initial compound).

DBP analysis

4 THMs, 4 haloacetonitriles (HANs), 2 haloketones (HKs) and chloropicrin were extracted following EPA method 551, which consists in a liquid-liquid extraction using MTBE, followed by detection using an Agilent 7890A gas chromatograph coupled with electron capture detector (GC-ECD) (Munch et al., 1995a). 9 HAAs were extracted and analyzed following EPA method 552.2, which is based on a liquid-liquid extraction with MTBE in acidic conditions followed by derivatization to methyl esters using acidic methanol (Munch et al., 1995b). Extracts were analysed using an Agilent 7890A GC system coupled with a 5975C mass spectrometer (GC-MS). HAcAms were analysed following the same protocol than EPA method 551, but using ethyl acetate for the liquid-liquid extraction step. All DBPs analysed are described in Table 2.

RESULTS AND DISCUSSION

Formation of haloacetamides from amino acids

The formation of haloacetamides was investigated after 72h of chlorination and chloramination (15 mg/L as Cl2) of 4 amino acids (50 µM each) at pH 8 (Figure 1). Asparagine was found to be a major dichloracetamide (DCAcAm) precursor during chloramination, producing up to 1160 µg/L DCAcAm. This can be attributed to the amide group present in the asparagine structure (Table 1). However, this does not explain why chloramination formed much more DCAcAm than chlorination. Other amino acids were much lower precursors of DCAcAm. In the case of tyrosine and tryptophan, more DCAcAm was formed by chlorination. These results indicate that different reaction pathways occur during chlorination and chloramination of amino acids. Especially, a specific pathway seems to occur during the chloramination of asparagine. [START_REF] Huang | Dichloroacetonitrile and dichloroacetamide can form independently during chlorination and chloramination of drinking waters, model organic matters, and wastewater effluents[END_REF] proposed a formation mechanism to explain the formation of DCAcAm from asparagine, but did not describe a specific pathway for NH2Cl, because the formation of DCAcAm was similar by chlorination and by chloramination in their experimental conditions (100 µM asparagine, 10 mg/L as Cl2, t = 2h, pH 6.9). were studied, each consisting in a different set of nutrients conditions (1: nitrates and phosphates added every day; 2: nitrates and phosphates added once; 3: nitrates, phosphates and silica added every day; 4: nitrates, phosphates and silica added once). All four samples were diluted with synthetic seawater to equal concentrations of DOC (1.5 mg C/L) and were chlorinated and chloraminated (5 mg Cl2/L) in the lab during 72h. Results showed a higher formation of haloacetamides by chlorination than by chloramination (Figure 3.a). The first mesocosm formed the highest concentration of HAcAms (240 µg/L), probably because of the presence of a higher concentration in dissolved organic nitrogen (DON = 0.19 mg N/L in samples from mesocosm 1 and 3; DON = 0.12 mg N/L in samples from mesocosms 2 and 4) due to continuous addition of nitrates. These results indicate that the nitrogen atom of haloacetamides would predominantly originate from organic nitrogen contained in algal cells, and not from NH2Cl. This implies that the predominant mechanism of HAcAms formation from algae cells is the decarboxylation pathway, where the nitrogen atom of HAcAms originates from DON. This is also suggested by the similar trend in HANs formation, exhibiting higher concentrations from chlorination than from chloramination (Figure 3.b). In the decarboxylation pathway, HAcAms are produced through HANs hydrolysis. Hence our results suggest a similar pathway for HANs and HAcAms formation, favoured by the presence of free chlorine. As it could be expected because of the presence of bromide ion in the seawater samples (60-65 mg/L as Br -), the proportion of brominated DBPs formed after chlorination and chloramination was more important than the amount of chlorinated DBPs. TBM, DBAcAm, DBAA, and DBAN were the predominant species among THMs, HAcAms, HAAs and HANs, respectively. This is explained by the fast reaction between chlorine and bromide ion leading to the production of hypobromous acid (HOBr) (Equation 2), which is an oxidizing species even more reactive than hypochlorous acid (HOCl) [START_REF] Grasshoff | Formation and distribution of haloacetic acids, THM and TOX in chlorination of bromide-rich lake water[END_REF].

HOCl + Br-→ HOBr + Cl-(2)

CONCLUSIONS

The aim of our work was to study the formation of HAcAms from chlorination and chloramination of marine algal organic matter. DBP formation potential tests were performed from model compounds (amino acids) and from simulated algal blooms in real seawater conditions. Results show that chlorination produces more HAcAms than chloramination from real algae (Synecococcus sp.), thus indicating that the nitrogen of HAcAms comes predominantly from DON through the decarboxylation of amino acids rather than from NH2Cl. Results of HAcAm formation potentials

Figure 1 .

 1 Figure 1. DCAcAm formation potential from chlorination and chloramination (15 mg Cl2/L) of amino acids (50 µM) after 72h contact time at pH8 (10 mM carbonate buffer).

Figure 2 .

 2 Figure 2. Formation of DCAcAm and TCAcAm from a) chlorination and b) chloramination of amino acids ([amino acids] = 50 µM; [oxidant] = 15 mg/L as Cl2; t = 72h; pH = 8).

Figure 3 .

 3 Figure 3. Formation of a) total HAcAms and b) total HANs from chloramination and chlorination (5 mg/L as Cl2) of four simulated algal blooms (mesocosms) (1.5 mg C/L) during 72h.

Table 2

 2

	.Description of monitored DBPs
	Method of	Class of compounds and	Abbreviation	Full name
	determination	structures
	EPA 551	Trihalomethanes (THMs)	TCM	Trichloromethane (chloroform)
			DCBM	Dichlorobromomethane
			DBCM	Dibromochloromethane
			TBM	Tribromomethane (bromoform)
		Haloacetonitriles (HANs)	DCAN	Dichloroacetonitrile
			TCAN	Trichloroacetonitrile
			BCAN	Bromochloroacetonitrile
			DBAN	Dibromoacetonitrile
		Haloketones (HKs)	1,1-DCP	1,1-dichloropropanone
			1,1,1-TCP	1,1,1-trichloropropanone
		Halonitromethane	TCNM	Trichloronitromethane
	EPA 551 -EA Haloacetamides (HAcAms)	CAcAm	Chloroacetamide
			BAcAm	Bromoacetamide
			DCAcAm	Dichloroacetamide
			TCAcAm	Trichloroacetamide
			DBAcAm	Dibromoacetamide
			TBAcAm	Tribromoacetamide
			BCAcAm	Bromochloroacetamide
			CIAcAm	Chloroiodoacetamide
			BIAcAm	Bromoiodoacetamide
			DIAcAm	Diiodoacetamide
	EPA 552.2	Haloacetic acids (HAAs)	MCAA	Monochloroacetic acid
			MBAA	Monobromoacetic acid
			DCAA	Dichloroacetic acid
			TCAA	Trichloroacetic acid
			BCAA	Bromochloroacetic acid
			DCBAA	Dichlorobromoacetic acid
			DBAA	Dibromoacetic acid
			DBCAA	Dibromochloroacetic acid
			TBAA	Tribromoacetic acid

from amino acids are in agreement with this pathway, tyrosine and tryptophan exhibiting higher DCAcAm formation by chlorination than by chloramination. However, asparagine was found to be a major DCAcAm precursor during chloramination, thus indicating a specific reaction pathway. While the amide group of asparagine is probably the main reason of its important DCAcAm formation potential, the specific role of NH2Cl as compared to HOCl needs to be further investigated. The use of isotopically labelled monochloramine( 15 NH2Cl) could help understanding the incorporation of nitrogen from asparagine or from NH2Cl. Dihalogenated species of DBPs (i.e., DBAcAm, DBAA, DBAN) were the dominant species formed by both chlorination and chloramination of mesocosms samples, while chlorination of amino acids formed more TCAcAm than DCAcAm, indicating a higher incorporation of chlorine.

Future work will focus on the formation of haloacetamides by chloramination of different marine AlOM obtained from algae monocultures grown in laboratory to understand the specificities of haloacetamides precursors from real algae. Reactions will be conducted in the presence isotopically labelled 15 NH2Cl to assess the proportion of nitrogen coming from organic nitrogen sources and from monochloramine.