Embedding Anatomical or Functional Knowledge in Whole-Brain Multiple Kernel Learning Models

Pattern recognition models have been increasingly applied to neuroimaging data over the last two decades. These applications have ranged from cognitive neuroscience to clinical problems. A common limitation of these approaches is that they do not incorporate previous knowledge about the brain structure and function into the models. Previous knowledge can be embedded into pattern recognition models by imposing a grouping structure based on anatomically or functionally defined brain regions. In this work, we present a novel approach that uses group sparsity to model the whole brain multivariate pattern as a combination of regional patterns. More specifically, we use a sparse version of Multiple Kernel Learning (MKL) to simultaneously learn the contribution of each brain region, previously defined by an atlas, to the decision function. Our application of MKL provides two beneficial features: (1) it can lead to improved overall generalisation performance when the grouping structure imposed by the atlas is consistent with the data; (2) it can identify a subset of relevant brain regions for the predictive model. In order to investigate the effect of the grouping in the proposed MKL approach we compared the results of three different atlases using three different datasets. The method has been implemented in the new version of the open-source Pattern Recognition for Neuroimaging Toolbox (PRoNTo).


Introduction
During the last years there has been a substantial increase in the application of machine learning models to analyse neuroimaging data (please see Pereira et al. 2009 andHaynes 2015 for overviews). In cognitive neuroscience, applications of these models -also known as brain decoding or mind reading-aim at associating a particular cognitive, behavioural or perceptual state to specific patterns of brain activity. In the context of clinical neuroscience, machine learning analyses usually focus on predicting a group membership (e.g. patients vs. healthy subjects) from patterns of brain activation/anatomy over a set of voxels. Due to their multivariate properties, these approaches can achieve relatively greater sensitivity and are therefore able to detect subtle and spatially distributed effects. Recent applications of machine learning models to neuroimaging data include predicting, from individual brain activity, the patterns of perceived objects (Haynes and Rees 2005;Ramirez et al. 2014), mental states related to memory retrieval (Polyn et al. 2005) and consolidation (Tambini and Davachi 2013), hidden intentions (Haynes et al. 2007) and semiconstrained brain activity (Schrouff et al. 2012). These techniques also showed promising results in clinical applications J. Schrouff and J. M. Monteiro equally contributed to the present work.
* Jessica Schrouff jschrouff@ucl.ac.uk 1 (see e.g. Klöppel et al. 2012), providing potential means of computer-aided diagnostic tools for Alzheimer's disease (Klöppel et al. 2008), Parkinson's disease (e.g. Orrù et al. 2012;Garraux et al. 2013) or depression (Fu et al. 2008). Accordingly, various software packages have been implemented to ease the application of machine learning techniques to neuroimaging data. To cite a few: The Decoding Toolbox (Hebart et al. 2015), MVPA toolbox, PyMVPA (Hanke et al. 2009a, b), Nilearn (Abraham et al. 2014), Representational Similarity Analysis (Kriegeskorte et al. 2008), CoSMoMVPA (Oosterhof et al. 2016), Searchmight Botvinick 2011), 3Dsvm (LaConte et al. 2005), Probid, Mania (Grotegerd et al. 2014), PETRA or our own work PRoNTo (Schrouff et al. 2013a). When applying machine learning predictive models to whole brain neuroimaging data a researcher often wants to be able to answer two questions: (1) Which brain regions are informative for the prediction? (2) Why are these regions informative? Considering the first question, although linear models generate weights for each voxel, the model predictions are based on the whole pattern and therefore one cannot arbitrarily threshold the weights to identify a set of informative features (or voxels). Indeed, if one were to threshold a weight map (e.g. by removing voxels/regions with low contribution), the result would be a new predictive function that has not been evaluated. In order to identify which features have predictive information one can use feature selection approaches or sparse models. One limitation of these approaches is that often they do not take into account our previous knowledge about the brain. We know that the brain is organised in regions and the signal within these regions are expected to vary smoothly. One way to incorporate this knowledge into the models is to use structured or group sparsity. A number of studies have shown the benefits of using structured sparse approaches in neuroimaging applications (e.g. Baldassarre et al. 2012, Grosenick et al. 2011. However, these models are computationally expensive and it is difficult to design a structured sparsity that incorporates all characteristics of the neuroimaging data. An alternative way to incorporate knowledge about the data into the models is to use group sparsity. For example, there is evidence that group sparse regularization (i.e. group lasso) can improve recovery of the model's coefficients/weights in comparison with the lasso when the grouping structure is consistent with the data (Huang and Zhang 2010). Here, we used anatomical/functional information to define the grouping structure and a sparse version of Multiple Kernel Learning (MKL) to simultaneously learn the contribution of each brain region to the predictive model.
The question of why a set of regions carries predictive information is more difficult to answer and has been previously discussed in the literature (e.g. Haufe et al. 2014;Weichwald et al. 2015;Kia et al. 2016). Basically, weights of linear predictive models show the relative contribution of the features for prediction, but do not disentangle potential causes for the contribution. For example, as shown by Haufe and collaborators (Haufe et al. 2014), a feature might have a high weight (or a high contribution) due to an association with the labels or a high weight to cancel correlated noise between the features. Therefore, we argue that additional analysis needs to be done (e.g. univariate statistical tests) to understand why a specific feature (or region) has a high contribution to a predictive model.
In this work, we propose an approach that is able to select a subset of informative regions for prediction based on an anatomical/functional atlas, thereby addressing the first question. However, we do not attempt to address the second question, as we believe multivariate predictive models cannot provide a clear answer to why a specific feature/region has a high contribution to the model (Weichwald et al. 2015). In the present work, we will refer to the 'interpretability' of a predictive model as its ability to identify a subset of informative features/regions.

Related Approaches
Different solutions have been proposed to identify which features contribute to the model's prediction 1 : Kriegeskorte et al. (2006) proposed a locally multivariate approach, known as Bsearchlight^, whereby only one voxel and its direct neighbours (within a sphere which radius is defined a priori) are selected to build the machine learning model. This operation is then repeated for all voxels, leading to a map of performance (e.g. accuracy for classification and mean squared error, MSE, for regression). Based on the significance of model performance in each sphere, the resulting maps can be thresholded. While this approach can provide insights on which regions in the brain have a local informative pattern, it presents the disadvantage of considering each sphere independently. The brain is therefore not considered as a whole anymore, but as a collection of partially overlapping spheres, which reduces the multivariate power of machine learning models by focusing only on local patterns. The interested reader can refer to Etzel et al. 2013 for a discussion on the promise, pitfalls and potential of this technique.
Another approach that has been used to threshold weight maps is to perform a permutation test at each voxel to generate a map of p-values (e.g. Mourão-Miranda et al. 2005, Klöppel et al. 2008. In this case the labels of the training data are randomly shuffled p times and the model is trained using the shuffled labels to generate a null distribution of the models' weight for each voxel. The voxels with a statistically high contribution (positive or negative) to the model compared to its null distribution can then be highlighted. The resulting statistical maps can be thresholded, using the p-values obtained for each voxel. The correction for multiple comparisons should be performed with care, as detailed in (Gaonkar and Davatzikos 2012). In addition, this approach is computationally expensive.
Some authors have proposed the use of sparse models, like LASSO (Tibshirani 1996) or Elastic-net (Zou and Hastie 2005), as they are able to estimate solutions for which only few voxels are considered relevant. Structured sparse models, such as sparse Total Variation (TV, Baldassarre et al. 2012) and Graph Laplacian Elastic Net (GraphNET, Grosenick et al. 2011), allow incorporation of domain knowledge through additional spatial and temporal constraints and carry the promise of being more interpretable than non-structured sparse methods, such as LASSO or Elastic Net methods. A drawback of the sparse models is that the solution is highly dependent on the way the prior or regularization term is specified. Often models with different regularization terms (e.g. LASSO, Elastic-net, Total Variation) achieve similar accuracies for different solutions (Baldassarre et al. 2012). In this sense, some authors have argued that the quality of spatial patterns extracted from sparse models cannot be assessed purely by focusing on prediction accuracy (Rasmussen et al. 2012).
Feature selection based on stability theory (Meinshausen and Bühlmann 2010) has also been proposed as a mapping approach by identifying a subset of stable features that are relevant to the predictive model (Rondina et al. 2014). This approach relies on the idea of choosing relevant features that are stable under data perturbation. Data are perturbed by iteratively sub-sampling both features and examples. For each perturbation, a sparse method (e.g. LASSO) is applied to a sub-sample of the data. After a large number of iterations, all features that were selected in a large fraction of the perturbations are selected. Although this approach has the potential to identify reliable relevant features for the predictive models, it does not account for prior knowledge about brain anatomy neither for the spatial correlation among the voxels.
Another approach to tackle the interpretability of machine learning models is to use previous knowledge about brain anatomy to segment the whole brain multivariate pattern into regional patterns. This strategy was used in (Schrouff et al. 2013b): the authors proposed local averages of the model weights according to regions defined by the Automated Anatomical Labelling (AAL, Tzourio-Mazoyer et al. 2002) atlas. Regions were then sorted according to their proportional contribution to the weight vector or decision function, thereby providing a ranking of the regions. Even though the results of this study showed that regions ranked in the top 10 (arbitrarily fixed threshold) were in line with previous univariate studies, this approach does not solve the issue of thresholding since for non-sparse machine learning models 2 (e.g. Support Vector Machines, Kernel Ridge Regression) all brain regions considered will have some contribution to the model's predictions. Investigating regional contribution through post-hoc summarization was also performed in Hanke et al. 2009a. In their work, the authors matched probabilistic weight maps with anatomical information to derive a 'specificity' measure for each region of interest. This approach however suffers from the same limitation, i.e. regions with low sensitivity are part of the decision function and cannot be pruned.
Multiple Kernel Learning (MKL, Bach et al. 2004) approaches have been previously applied in the context of neuroimaging to e.g. perform multi-modal diagnosis of Alzheimer disorders (Hinrichs et al. 2011;Zhang et al. 2011), attention deficit hyperactivity disorder (ADHD) children (Dai et al. 2012), predict cognitive decline in older adults (Filipovych et al. 2011) and discriminate three Parkinsonian neurological disorders (Filippone et al. 2012). In (Filippone et al. 2012), each kernel corresponded to either an image modality or an anatomically labelled region. The authors used the kernel weights to analyze the relative informativeness of different image modalities and brain regions. However, the considered algorithm was not sparse in the kernel combination, making it difficult to determine a subset of regions with highest contribution to the model. Our work differs from Filippone et al. 2012 as we use MKL as an exploratory approach to find a (sparse) subset of informative regions for a predictive model, considering all brain regions a priori defined by a whole brain template.

Proposed Approach
The proposed framework combines anatomical/functional parcellations of the brain, MKL and sparsity. More specifically, we use a sparse version of the MKL algorithm to simultaneously learn the contribution of each brain region, previously defined by an atlas, to the decision function. As the considered technique is sparse, some kernels (here corresponding to brain regions) will have a perfectly null contribution to the final decision function. The resulting weight maps at the voxel and region levels will hence be sparse and do not need to be thresholded. In summary, here we investigate the introduction of anatomical or functional a priori knowledge in a MKL whole brain model and compare the results when using different atlases, both in terms of model performance and obtained weight maps. The proposed approach has two potential benefits: (1) it can lead to improved overall generalisation performance when the grouping structure imposed by the atlas is consistent with the data; (2) it can identify a subset of relevant brain regions for the predictive model. It is important to note that our approach does not provide information about why a specific feature has a high weight (or contribution) to the model (Haufe et al. 2014, Weichwald et al. 2015 but rather aims at identifying a (sparse) list of regions that contribute to the model's predictive function. The approach is illustrated using three different atlases (described in the methods section) and three public datasets: the functional MRI (fMRI) Haxby dataset (Haxby et al. 2001) which investigates the differences in brain activity when viewing different types of visual stimuli, the fMRI 'face' data set (Henson et al. 2002) which studies changes in brain activity when looking at images of faces (famous, nonfamous and scrambled), and the structural MRI (sMRI) OASIS dataset (Open-Access Series of Studies, oasis-brains.org; Marcus et al. 2007), which consists of structural images obtained from non-demented and demented older adults. The method was implemented in PRoNTo (http://www.mlnl.cs.ucl.ac.uk/pronto/).

Datasets and Pre-Processing
Three public datasets were used to illustrate the proposed approach. The first one has been previously used in pattern recognition for neuroimaging studies (Haxby et al. 2001;Hanson et al. 2004;O'Toole et al. 2005) and for describing the functionalities of different software toolboxes (Hanke et al. 2009a, b;Schrouff et al. 2013a). The data consist of a block design fMRI experiment acquired using a visual paradigm, where the participants passively viewed grey scale images of eight categories: pictures of faces, cats, houses, chairs, scissors, shoes, bottles, and control, non-sense images. As an illustrative example, we chose to analyse the data from a single subject (participant 1), consisting of 12 runs, each comprising eight blocks of 24 s showing one of the eight different object types and separated by periods of rest. Each image was shown for 500 ms followed by a 1500 ms inter-stimulus interval. Fullbrain fMRI data were recorded with a volume repetition time of 2.5 s. Each category block therefore corresponds approximately to nine scans, separated by six scans of rest. For further information on the acquisition parameters, please consult the original reference (Haxby et al. 2001). The data were preprocessed using SPM8 (http://www. Fil.ion.ucl.ac.uk/spm/ software/). We motion corrected, segmented and normalized the scans according to the MNI template. No smoothing was applied to the data. 3 For proof of concept, we chose to focus the analysis on the comparison between viewing 'faces' and viewing 'houses', since it was reported as leading to high accuracy values and precise anatomical localization of the most discriminative regions (Schrouff et al. 2013a). Therefore we expected visual areas to have a high contribution to the predictive model.
The second dataset consisted of a single subject eventrelated fMRI data freely available from the SPM website comprising a repetition priming experiment, where two sets of 26 familiar (famous) and unfamiliar (non-famous) faces were presented against a checkerboard baseline. A random sequence of two presentations of each face was created from each set. The faces were presented for 500 ms with a stochastic distribution of stimulus onset asynchrony (SOA) determined by a minimal SOA of 4.5 s and 52 randomly interspersed null events. The subject was asked to make fame judgments by making key presses. Whole brain fMRI data were recorded with a volume repetition time of 2 s. For further information on the acquisition parameters, please consult the original work (Henson et al. 2002). The data were pre-processed using SPM8. This included motion correction, segmentation, normalization to the MNI template and smoothing ([8 8 8] mm). To classify famous versus non-famous faces we first fitted a GLM to all voxels within the brain, using SPM8. The design matrix comprised as many columns as events (all famous and non-famous faces presented, in order to obtain one beta image per event) plus the movement parameters and the mean regressor. The betas corresponding to the second repetition of famous and non-famous faces were used for classification.
The third dataset, the Open-Access Series of Studies (OASIS, oasis-brains.org; Marcus et al. 2007), illustrates the potential of the proposed methodologies in clinical settings. It consists of structural MRI images from nondemented and demented older adults. In the OASIS dataset, patients were diagnosed with dementia using the Clinical Dementia Rating (CDR) scale as either nondemented or with very mild to mild Alzheimer's disease (Morris 1993). A global CDR of 0 indicates no dementia (healthy subjects) and a CDR of 0.5, 1, 2 and 3 represent very mild, mild, moderate and severe dementia respectively. The patients were age and gender matched with the controls, such that our analysis comprises the structural MRI images from fifty patients diagnosed with very mild and mild dementia (M = 75.3, SD = 6.8, 28 females) and fifty healthy controls (M = 75, SD = 6.7, 28 females). The OASIS data were also pre-processed using SPM8. The first step was to average all the repeats for each session followed by a grey matter segmentation, then, the segmented images were normalized and smoothed with a Gaussian kernel with a full width at half maximum (FWHM) of [8 8 8] mm.
Additional pre-processing was applied before the machine learning modelling. The data were linearly detrended (fMRI data only, polynomial detrend of order 1). In order to ensure that the MKL and SVM models were based on the same set of voxels we built one binary mask from each atlas. In the case of the fMRI data, the mask defined by the considered atlas, was applied to each image to select the voxels used as a feature in the modelling. In the case of the structural MRI data, we first selected voxels that had a probability of being located in grey matter equal or above 30% in all subjects and then applied a mask defined by the considered atlas to select the voxels. For all datasets a linear kernel was built for each region as defined by the considered atlas.
Three atlases were used to investigate the effect of using different anatomical or functional priors in the MKL model ( Fig. 1): (1) The Automated Anatomical Labeling (AAL, Tzourio-Mazoyer et al., 2002) atlas, built using the WFU-PickUp Atlas toolbox of SPM and consisting of 116 brain regions. This atlas is a widely used manual macroanatomical parcellation of the single subject MNI-space template brain.
(2) The Brodmann + atlas, built using the WFU-PickUp Atlas toolbox of SPM and consisting of 75 regions. This atlas includes 47 out of the 52 areas defined by K. Brodmann, based on cytoarchitecture or histological structure, as well as other structures and nuclei.
(3) The atlas built from the Human Connectome Project (HCP, Glasser et al. 2016). This multi-modal parcellation (atlas) is probably the most detailed cortical in-vivo parcellation available to date. The HCP MMP 1.0 has been built using surface-based registrations of multimodal MR acquisitions and an objective semi-automated neuroanatomical approach to delineate 180 areas per hemisphere bounded by sharp changes in cortical architecture, function, connectivity, and/or topography in a group average of 210 healthy young adults from the HCP cohort. It comprises 180 bilateral regions.
In all MKL models, the kernels were mean centred and normalized before classification, taking the training set/test set split into account. Mean centring the kernel corresponds to mean centre the features across samples (i.e. it is equivalent to subtracting the mean of each feature/voxel, computing the mean based on the training data), while normalizing the kernel corresponds to dividing each feature vector (i.e. each sample) by its norm. The later operation is particularly important when using MKL approaches to compensate for the fact that the different kernels might be computed from different numbers of features (i.e. different region sizes). This operation can hence be seen as giving an equal chance to all regions, independently of their sizes. Both operations were considered as preprocessing steps and can affect the model performance and obtained weight maps. For single kernel modelling (i.e. SVM models), the kernels were first added to provide the whole brain feature set, then mean centred. It should be noted that adding linear kernels is equivalent to concatenating the features/voxels. The resulting kernel is not normalized as this is not a common operation for single kernel modeling and often leads to decreases in model performance (unpublished results).

Machine Learning, Modelling
The two classifiers considered in the present work are based on binary SVM machines (Boser et al., 1992). More specifically, single kernel analyses were conducted using the LIBSVM implementation of SVM (Chang and Lin 2011), while multi-kernel learning was performed using the SimpleMKL package (Rakotomamonjy et al., 2008), which resorts to the SimpleSVM algorithm (Canu et al., 2003). The framework of those two procedures is described below:

Single Kernel Modelling
Mathematically, let X∈R n,l , the data matrix of samples (n) by features (l) and y∈R n the corresponding labels, where each row x i corresponds to a feature vector and y i corresponds to its respective label. Supervised learning approaches for binary classification, such as the SVM, estimate a decision function f, which separates the data into different classes defined by the labels. In a linear model, f is of the form (Eq. 2.1): With <,> representing the dot product between the weight vector w ∈ R l and a feature vector x i , and b being a bias term.
The decision function f of an SVM is obtained by solving the following optimisation problem (Boser et al., 1992): Where i indexes the samples, from 1 to n, C corresponds to the soft-margin parameter, ∑ i ξ i is an upper-bound on the number of training errors and b is a bias term. The solution of the optimisation problem can be written as (please see appendix A1 for details): Substituting w into Eq. 2.1 and considering the linear kernel definition K(x, x i ) = x, x i , we can re-write the decision function in its dual form as Where α i and b represent the coefficients to be learned from the examples and K(x, x i ), the kernel, is a function characterising the similarity between samples x and x i .
An illustration of whole brain single kernel modelling is presented in Fig. 2.

Multiple Kernel Learning
In multiple kernel learning, the kernel K(x, x ′ ) can be considered as a linear combination of M Bbasis^kernels , i.e.: Therefore, the decision function of an MKL problem can be expressed in the form: The considered multiple kernel learning approach is based on the primal formulation of an SVM binary classifier (Rakotomamonjy et al., 2008) and the solution can be obtained by solving the following optimisation problem: The computed kernels K m are added to obtain a whole brain linear kernel K. The kernel and its associated labels are used to train the model and estimate the model parameters w. The model can then be applied to new/unseen data x* to obtain an associated predicted label With d m representing the contribution of each kernel K m to the model. Therefore, both d m and w m have to be learned simultaneously. In this formulation, proposed by (Rakotomamonjy et al., 2008), the L1 constraint on d m enforces sparsity on the kernels with a contribution to the model. Furthermore, it results in a convex optimisation problem that can be solved using a simple SVM machine on K and gradient descents to find d m . For further details, please refer to (Rakotomamonjy, et al., 2008).
For the considered MKL optimisation problem the weights w m can be expressed as (please see appendix A1 for details) In the present case, MKL can be seen as a feature selection technique, i.e. each kernel corresponds to a different subset of features (corresponding to the labelled regions). The considered approach is illustrated in Fig. 3. However, MKL can potentially be used as a model selection strategy, where each kernel corresponds to a different model (e.g. different parameter of a non-linear kernel, Rakotomamonjy et al., 2008). In a neuroimaging context, MKL approaches were mostly used to combine heterogeneous sources of features, such as different imaging modalities (e.g. Filippone et al., 2012) or imaging with psychological testing (e.g. Filipovych et al., 2011). Such combination of multiple image modalities can also be performed using the MKL implementation in PRoNTo v2.0.

Assessing Performance
We performed a nested cross-validation procedure to train the model and optimise the model's hyperparameters. The external loop was used for assessing the model's performance and the internal loop was used for optimising the models hyperparameters (soft-margin parameter, C, for the SVM and SimpleMKL). For all models (MKL and SVM) the hyperparameter range was [0.01, 1, 100]. The reason for the limited number of tested values was the high computational cost of MKL with parameter optimisation. For the Haxby dataset we used a leave-one-block-out cross-validation for the external loop and the internal loop. For the 'face' dataset, we performed a leave-one-example-per-class-out cross-validation, for the external and internal loop. For the OASIS dataset we used a k-folds cross-validation on subjects-pergroup, with k = 10 folds for the external loop (i.e. leaving 10% of the subjects out, half of them being demented, half being healthy) and k = 5 folds for the internal loop. Model . The model can then be applied to new/unseen data x* to obtain an associated predicted label, based on feature vectors defined using the same atlas, x 1 *, x 2 *, …, x M * performance was assessed by balanced accuracy values, computed as the average of the class accuracies (corresponding to the sensitivity and specificity). A p-value was associated to each accuracy measure using permutation tests: the labels of the examples in the training set were randomly shuffled (taking the block structure of the datasets into account) before building a model. Results were considered significant when the obtained models performed equally or better than the model without shuffling the labels at most 5% of the time across 100 permutations (i.e. p-value < 0.05).

Weight Map
As shown in Eq. 2.8 the models weights (w and w m ), representing the contribution of each feature (here voxel) for the decision function or predictive model can be explicitly computed and plotted as brain images in order to display the decision function of the model based on previously defined brain regions. To avoid scaling issues between weight maps (e.g. from different folds or data sets), the resulting weight maps were normalized (i.e. w/||w|| 2 ).
As our MKL approach can be seen as a hierarchical model of the brain, it is possible to derive weights at two levels: (1) the weight or contribution of each region to the decision function, i.e. the values of d m , and (2) the weights for each voxel (see appendix A2 for the derivation of the weights per voxel). The weights at the voxel level can provide insights on the homogeneity of the discriminative patterns within the regions. Regions were ranked according to their contribution to the model (i.e. d m ), averaged across folds. Only regions with a positive (i.e. non-null) contribution to the decision function f are displayed (i.e. #d m > 0).

Stability of the Regions' Contribution
To investigate whether the selected regions are stable across the folds of the cross-validation (i.e. variability in the training data), we computed the Breproducibility^of the regions' ranking. Firstly, the ranking of a region is computed within each fold by sorting the kernel contributions in ascending order. Regions with a null contribution were assigned a null rank. The minimum value of the ranking is hence 0, while its maximum corresponds to the number of regions. The Expected Ranking (ER) is computed as the average of the ranking across folds. As in (Kia et al. 2016), we compute the cosine of the angle between the expected ranking (ER) and the ranking in each fold and estimate the 'reproducibility' of the ranking as the expectation of the cosine. This measure provides an estimation of the 'distance' between the ranking in each fold and the average ranking.
More specifically, if we assume an angle α j between ER and R j , the ranking in fold j (j = 1… number of folds), we have (Eq. 2.9): The reproducibility ψ R of the ranking (0 < =ψ R < =1) is then (Eq. 2.10): The closer this number is to 1, the more stable the solution is across folds. It is important to note that these values are meaningful only if the corresponding model performs significantly above chance level.

Comparison of Atlases
We finally compare different priors (i.e. atlases) in terms of obtained weight maps. To this end, we computed the Pearson correlation between the weight maps at the voxel level of each atlas, for overlapping voxels (i.e. voxels considered for modeling in both atlases). We then obtained three values of correlation, one for each pair of atlas. The closer this value is to one, the more similar the two considered weight maps are. The significance of the obtained correlation values was tested using 1000 non-parametric permutations. As correlation measures do not take into account null values, we also estimated the proportion of null weights that is shared by both atlases (i.e. voxels with 0 weight in both atlases, the intersection of null values) compared to the total number of overlapping voxels.

Haxby Dataset
Model Performance Table 1 shows that the model can discriminate with high accuracy if the subject was viewing images of faces versus images of buildings, for all models and atlases. This was expected in view of the previous performances obtained using this dataset (e.g. Hanke et al. 2009a;Schrouff et al. 2013a). Overall, the MKL models perform better than the SVM models, with the MKL-HCP model leading to the best performance. For both MKL and SVM, the Brodmann atlas leads to a slight decrease in balanced accuracy when compared to the AAL and HCP atlases. Please note that using only the left and right fusiform regions as defined by the AAL atlas leads to a balanced accuracy of 99.5% (108/108, 107/108). This shows that these visual areas carry a lot of predictive information. These regions are therefore expected to have a high model contribution and ER.

Stability of the Regions' Contribution
For each MKL model, we present the number of regions selected (i.e. with a non-null contribution across folds) in Table 2, as well as the model's reproducibility.
For this dataset, all models are quite sparse, with a relatively low number of regions with a non-null contribution to the model. The models with the highest accuracies (namely AAL and HCP) also lead to the highest reproducibility.

Comparison of Weight Maps Across Atlases
The weight maps for each atlas (at the voxel level) are displayed in Fig. 4 and the list of selected regions with non-null contributions for the MKL models for each atlas are displayed in appendices Tables 7, 8 and 9, along with their contributions d m and expected ranking ER. We can see that the fusiform regions (left and right) are ranked highly in the MKL-AAL model (ranks 115/116 and 103/ 116, respectively). Similarly, the MKL-Brodmann model selected area 19 (visual cortex, V3, V4 and V5) with highest rank (70/74), and area 37 (overlapping with the fusiform gyrus) with rank (50/74). In contrast, the MKL-HCP model selected ventromedial areas 1 and 2 with ranks (180/180) and (170/180), respectively.
In order to verify if the weight maps for the different MKL models were similar we computed the pairwise correlation coefficient between the weight vectors for the different models. The weight vectors for the AAL-MKL and Brodmann-MKL models have a correlation coefficient of ρ = 0.5330 (p = 9.9e −4 ), with 79.16% of voxels with a null weight in both models. The AAL-MKL and HCP-MKL weight vectors have a correlation coefficient of ρ = 0.3470 (p = 9.9e −4 ) and shared 83.05% of null weights, while the Brodmann-MKL and HCP-MKL weight vectors have a correlation coefficient of ρ = 0.5402 (p = 9.9e −4 ) with 86.51% of common null weights. The weight vectors for the two models leading to the highest performance and reproducibility are hence significantly correlated. Table 3 displays model performance for the MKL and SVM models considered. Most of the models were able to discriminate if the subjects were looking at 'famous' vs 'non-famous faces', however the accuracies were lower than the ones observed for the Haxby dataset. For this dataset there is an improvement in performance for the MKL models based on the AAL and Brodmann atlases with respect to the SVM models. The Brodmann-MKL model has the best performance across the MKL models. Results for the HCP-MKL and AAL-SVM models are not significant.

Stability of the Regions' Contribution
For each MKL model, we present the number of regions selected (i.e. with a non-null contribution across folds) in Table 4, as well as the model's reproducibility.
For this dataset, between 35% and 48% of the regions were selected, resulting in moderate sparsity. As for the Haxby dataset, the atlases leading to the best performance (namely AAL and Brodmann) lead to the highest reproducibility.

Comparison of Weight Maps Across Atlases
The weight maps for each atlas (at the voxel level) are displayed in Fig. 5 and the list of selected regions with nonnull contributions for the MKL models for each atlas are displayed in appendices Tables 10, 11 and 12, along with their True positives (resp. negatives) represent the class accuracy for faces (resp. houses) samples classified correctly as faces (resp. houses). Note that the difference between the SVM models is only the mask used to select the voxels, which is based on the atlas The correlation coefficient between the weight vectors for the AAL-MKL and Brodmann-MKL models is ρ = 0.1717 (p = 9.9e −4 ), with 34.23% of voxels with a null weight in both atlases. The AAL-MKL and HCP-MKL weight vectors have a correlation coefficient of ρ = 0.2760 (p = 9.9e −4 ) and shared 44.43% of null weights, while the Brodmann-MKL and HCP-MKL weight vectors have a correlation coefficient of ρ = 0.1550 (p = 9.9e −4 ) with 33.00% of common null weights. For this dataset, the similarity between weight maps is much lower than for the Haxby dataset, with most null weights being so in only one atlas.

Model Performance
Classifying healthy versus demented patients (with mild and very mild dementia) led to the accuracy values presented in Table 5. All models led to significant classification results. The results show that SVM models perform better than the  True positives (resp. negatives) represent the class accuracy for 'famous faces' (resp. 'non-famous faces') samples classified correctly as 'famous faces' (resp. 'non-famous faces'). Note that the difference between the SVM models is only the mask used to select the voxels, which is based on the atlas MKL models for this dataset, and that the Brodmann atlas led to highest performance for both MKL and SVM models.

Stability of the Regions' Contribution
For each MKL model, we present the number of regions selected (i.e. with a non-null contribution across folds) in Table 6, as well as the model's reproducibility.
The decision function seems to be based on a more distributed set of regions for this dataset. This was also supported by the higher model performance of SVM compared to MKL, since the SVM is a non-sparse model. As observed in the other datasets, the model leading to the highest accuracy (i.e. using Brodmann atlas) leads to the highest reproducibility.

Comparison of Weight Maps Across Atlases
The weight maps for each atlas (at the voxel level) are displayed in Fig. 6 and the list of selected regions with non-null contributions for the MKL models for each atlas are displayed in appendices Tables 13, 14  The correlation coefficient between the weight vectors at the voxel level for the AAL-MKL and Brodmann-MKL models is ρ = 0.4213 (p = 9.9e −4 ), with 7.79% of voxels with null weights in both atlases. The AAL-MKL and HCP-MKL weight vectors have a correlation coefficient of ρ = 0.3507 (p = 9.9e −4 ) and shared 18.15% of null weights, while the Brodmann-MKL and HCP-weight vectors have a correlation coefficient of ρ = 0.4632 (p = 9.9e −4 ) with 9.43% of common null weights. For this dataset, there were more similarities between the AAL-MKL and the Brodmann-MKL models and between the Brodmann-MKL and the HCP-MKL models than between the AAL-MKL and the HCP-MKL models, which both have lower accuracies.

Discussion
In this work, we present a novel approach to introduce anatomical or functional information in whole-brain machine learning models. Our procedure combines a priori information about the brain anatomy or function from an atlas with Multiple Kernel Learning (MKL, Rakotomamonjy et al., 2008), thereby estimating the contribution of each previously defined region of interest for the predictive model. Furthermore, the considered algorithm is sparse in the number of kernels (L1-norm constraint), therefore it selects a subset of regions that carry predictive information. Our approach results in a list of pre-defined brain regions, which can be ranked according to their contribution to the model. As previously mentioned, the obtained list of regions does not need to be thresholded since the regions which were not selected by the model in any fold have a null contribution to the model (i.e. d m = 0). This is a clear asset over techniques such as summarising region weights post-hoc (Schrouff et al., 2013b) or the locally multivariate searchlight approach (Kriegeskorte et al., 2006). In the proposed approach, there is indeed no need to apply statistical tests to select regions with significant contributions and to apply corrections for multiple comparisons.
Our results show that the MKL model combining anatomically or functionally labelled regions had higher performance in comparison with the SVM model for the two fMRI datasets considered but not for the structural one (OASIS). These results suggest that the model assumptions of the MKL implementation considered (i.e. only a small subset of regions carry predictive information) is adequate for the fMRI datasets but not for the OASIS dataset. However, it is important to notice that the MKL results are also affected by the choice of the atlas (or the grouping structure), with some atlas being better than others depending on the dataset considered. The HCP atlas led to the best performance for the Haxby dataset, the AAL led to the best performance for the faces data set, while the Brodmann atlas led to best performance for the structural dataset. There was also a difference in performance for the SVM model as the voxels included for modeling are different from one atlas to another.
In terms of selected regions, the list of selected regions was dependent on the choice of the atlas. The different atlases considered have different brain coverage and very different parcellation of regions, therefore differences in the selection of regions are expected. In addition, the sparsity constraint of the MKL model might also contribute for the difference between the selected regions across atlases, as regions with correlated information might not be selected as being relevant for the predictive model.
As previously mentioned in the introduction, and discussed in previous works (Haufe et al., 2014, Weichwald et al., 2015, multivariate predictive models cannot provide a clear answer to why a specific region/feature has a high contribution to the model. Alternative approaches (e.g. univariate tests, correlation analyses, …) should be used to investigate why a set of regions has predictive information. In this work, we assess whether the highestranking regions are 'meaningful' by referring to the literature on the cognitive neuroscience and clinical issues they tackle.
Considering the Haxby dataset, the MKL models were able to discriminate with high accuracy if the subject was viewing images of faces versus images of buildings regardless of the atlas used. All MKL models identified regions that comprise the core system for visual processing as informative (Haxby et al., 2000 andHaxby et al. 2001). Nonetheless, the HCP atlas led to the best performance, reaching 100% balanced accuracy. The visual brain region with the highest contribution for HCP-MKL model was the ventro-medial visual area on the ventral surface of each hemisphere. It should be noted that this region has not been previously well parcellated in either the human or the macaque, usually being either left unparcellated or parts of it being included in other areas (Glasser et al., 2016). Surprisingly, the HCP-MKL was not able to select the fusiform face area. The way the visual system is differently parcellated in True positives (resp. negatives) represent the number of demented (resp. non-demented) patients classified correctly as demented (resp. non-demented). Note that the difference between the SVM models is only the mask used to select the voxels, which is based on the atlas For the face dataset, both the AAL-MKL and the Brodmann-MKL models were able to discriminate if the subjects were looking at famous faces versus unfamiliar faces. Surprisingly, the results for the HCP-MKL were not statistically significant according to the permutation test. A possible explanation for these results is the very small sample size (only 26 images per class), which can lead to very high variance in the model's performance, particularly when the leave-one-out cross validation framework is used (Varoquaux et al., 2017). Overall, the AAL-MKL model was able to identify regions that play an important role in visual processing such as cuneus, occipital regions and lingual gyrus and the Brodmann-MKL model was able to identify visuospatial area (BA7) and somatosensory association cortex (Minnesbusch et al., 2009;Liu et al. 2014). The selected regions for both atlases also included regions that have been implicated in recollection of episodic memories such as the precuneus and the posterior cingulate, (Nielson et al., 2010). In addition, prefrontal regions including the dorsolateral and ventromedial prefrontal cortex for Brodmann and the orbitofrontal cortex for both atlases were also selected. These regions have been found to be important for processing famous faces (Nielson et al., 2010;Sergerie et al. 2005;Leveroni et al., 2000), and may relate to the search and retrieval of person identity semantic information. The Brodmann-MKL model led to better performance than the AAL-MKL model. Interestingly, the posterior cingulate  Regions ranked further (i.e. rank >14) have a perfectly null contribution to the model. ER stands for Expected Ranking. The region size is displayed in voxels. The weights were averaged across folds (except for computing ER). A region label ending in 'L' (resp. 'R') means left (resp. right) hemisphere region cortex presented the highest weight in the Brodmann-MKL model with a contribution of 21.92% for the predictive model. The posterior cingulate has been implicated in recollection of episodic memories (Henson et al. 1999;Maddock et al. 2001), and consistently reported as involved in the processing of famous faces (Leveroni et al., 2000;Nielson et al., 2010) therefore it might have an important role in accessing information about famous people. In contrast, the Brodmann-MKL model attributed high weights to motor regions (precentral and midcingulate) suggesting that different types of faces might prompt different patterns of motor responses, which might be related to differences in reaction time between the two tasks (subjects were asked to make fame judgements during the two conditions, 'famous' and 'non-famous faces'). It should be noted that there were also similarities between the regions selected by the AAL-MKL and Brodmann-MKL, for example, the posterior cingulate cortex was selected by the AAL-MKL and the motor regions were selected by Brodmann-MKL. In summary, the face dataset presents a high variability between the top regions selected by the different MKL models depending on the atlas used. One possible explanation for this variability is that the classification task considered is more complex than the one presented in the Haxby dataset and therefore might involve a large network of regions, which are differently parcellated in the different atlases. For the OASIS dataset, our results show that the MKL models were able to discriminate with moderate accuracy between anatomical brain scans of patients with mild and very mild dementia and brain scans of age/gender-matched healthy controls. These results are in agreement with previous Alzheimer's literature, which shows that pattern recognition methods applied to structural MRI can consistently discriminate between brain scans of patient and healthy controls (Arimura et al., 2008, Klöppel et al., 2008, Magnin et al., 2009, Duchesne et al., 2008, Vemuri et al., 2008, Gerardin et al., 2009, Nho et al., 2010, Oliveira et al., 2010, Farhan et al., 2014. As expected, the MKL model was able to identify, for all atlases, regions that comprise the core system for episodic memories (including temporal regions, hippocampus, posterior cingulate and precuneus) and parieto-frontal regions, which are also in agreement with existing literature (Zhang et al., 2015, Magnin et al., 2009. The Brodmann atlas led to the highest performance, suggesting that this atlas has a good coverage and segmentation of the most discriminative regions. In the Brodmann-MKL model, the most informative region was the BA 7, which comprises the superior parietal lobule and part of precuneus. These regions have been described as diagnostic markers of Alzheimer's disease (Karas et al., 2007, Quiroz et al. 2013. Surprisingly, the Brodmann-MKL did not select the hippocampus as one of the most relevant regions for the prediction. In contrast, the hippocampus has been selected as the region with the highest contribution to the predictive model by the HCP-MKL and the left hippocampus has been selected as the third region with the highest contribution by the AAL-MKL. One possible explanation for the Brodmann-MKL model not selecting the hippocampus as a highly informative region might be due to the limitation of the considered MKL model, which might not select two regions as being important if they have correlated information. The considered MKL approach comprises a sparsity constraint on the kernels. This reflects an assumption that only a few kernels (or regions) carry predictive information, which might not be suited for all datasets. Our results show that the solutions for the three datasets considered had different degrees of sparsity. This might reflect that the pattern of interest is sparse for the Haxby and faces datasets, but not for the OASIS dataset. To dampen this issue, other regularization constrains, less conservative than the L1 could be envisaged. Future MKL developments including a combination of L1 and L2 regularizations should address this limitation. Our results also show that there was not an optimal atlas for all considered datasets. The definition of the atlas could also be approached as an optimisation problem, with the best grouping of regions being learnt automatically from the data.
In addition to identifying a subset of regions that contribute to the predictive model we can also investigate whether the selected results are stable across folds, i.e. is the selected subset of regions similar for slightly varying training sets? To investigate the stability of the selection of regions across folds, we used the expected ranking (ER) to estimate the reproducibility as a metric of stability. Across datasets, it seems that there is an association between reproducibility and how easy the discrimination between the categories is. This is illustrated by higher reproducibility score, ψ R, for the Haxby dataset (average across atlases: 0.9181), lower for the 'face' dataset (average across atlases: 0.8523) and the lowest for the OASIS dataset (average across atlases: 0.8174). These results suggest that, when the sparse constraint is appropriate for the classification problem, the L1 MKL model leads to both high performance and high reproducibility (e.g. Haxby dataset). For the face dataset, there seems to be balance between performance and reproducibility. This last result is in agreement with Kia et al. 2016, who showed a similar effect during model optimisation across most subjects of the MEG recordings of the face data set. For the OASIS dataset, the model leading to the highest performance also led to the highest reproducibility. However, as mentioned above, the sparse prior does not seem appropriate for this dataset, which prevents us from drawing further conclusions. A potential improvement of our approach would be to include the reproducibility of the expected ranking of regions as a criterion for optimising the soft-margin hyperparameter (C), in addition to the generalization performance. The advantages of introducing reproducibility as an additional optimisation criterion has been also discussed in (Rosa et al. 2015;Kia et al. 2016 andBaldassarre et al. 2017).

Implementation
The proposed MKL framework modelling the whole brain multivariate pattern as a combination of regional patterns has been implemented in our open-source software PRoNTo v2.0. The MKL algorithm corresponds to the simpleMKL toolbox (Rakotomamonjy et al., 2008). Detailed experiments on the memory usage and computational expenses can be found in the reference (Rakotomamonjy et al., 2008). Regarding memory use, the software (here Matlab) needs to be able to load all the ROI kernels simultaneously. Therefore, the size of the kernel (i.e. number of examples x number of examples) will play a role, as will the number of kernels built. Regarding CPU time, (Rakotomamonjy et al., 2008) showed that both the number of examples and the number of kernels affected computational expenses. For all the models considered in this work, running permutations to assess model significance was computationally expensive and was performed on a cluster for efficiency.
PRoNTo v2.0 also includes the post-hoc summarization of weights, as detailed in (Schrouff et al., 2013b). Nevertheless this approach suffers from the various limitations in terms of interpretation that were discussed in this work (i.e. the obtained list of regions should not be thresholded and the weights reflect the decision function of the model, not the neural sources of the signal).
Finally, PRoNTo v2.0 is provided under the GNU/GPL license 'as is' with no warranty. Our team has done its best to provide a robust framework to perform machine learning modeling of neuroimaging data. Such an endeavor is however a continuous process and we thank our users for reporting bugs. Improvements and bug fixes will be implemented in future versions of the software (v2.1 and v3.0, in progress).
In conclusion, here we present a new tool for introducing anatomical or functional information in whole-brain machine learning models using sparse multiple kernel learning. When the grouping structure defined by the atlas is consistent with the data and the sparsity constraint is appropriate, the proposed approach can lead to an increase in model performance when compared to whole-brain models. Furthermore, the obtained list of regions contributing to the model can then be investigated in terms of cognitive or clinical neuroscience and reproducibility.

Support Vector Machine (SVM)
The Lagrangian of the SVM problem (Eq. 2.2.) is given by Setting to zero the derivatives of the Lagrangian with respect to w m we get The Lagrangian of the MKL problem (Eq. 2.7.) is given by Setting to zero the derivatives of the Lagrangian with respect to w m we get

Haxby Data
Regions Selected in MKL-AAL Table 7 displays the regions with non-null contributions to the MKL-AAL model, along with their contribution d m , size (in voxels) and expected ranking ER, across folds. Table 8 displays the regions with non-null contributions to the MKL-Brodmann model, along with their contribution d m , size (in voxels) and expected ranking ER, across folds. Regions ranked further (i.e. rank >21) have a perfectly null contribution to the model. ER stands for Expected Ranking. The region size is displayed in voxels. The weights were averaged across folds (except for computing ER)

Regions Selected in MKL-Brodmann
Regions Selected in MKL-HCP Table 9 displays the regions with non-null contributions to the MKL-HCP model, along with their contribution d m , size (in voxels) and expected ranking ER, across folds.

Face Data
Regions Selected in MKL-AAL Table 10 displays the regions with non-null contributions to the MKL-AAL model, along with their contribution d m , size (in voxels) and expected ranking ER, across folds. Table 11 displays the regions with non-null contributions to the MKL-Brodmann model, along with their contribution d m , size (in voxels) and expected ranking ER, across folds. Table 12 displays the regions with non-null contributions to the MKL-HCP model, along with their contribution d m , size (in voxels) and expected ranking ER, across folds.

OASIS Data
Regions Selected in MKL-AAL Table 13 displays the regions with non-null contributions to the MKL-AAL model, along with their contribution d m , size (in voxels) and expected ranking ER, across folds. Regions ranked further (i.e. rank >13) have a perfectly null contribution to the model. ER stands for Expected Ranking. The region size is displayed in voxels. The weights were averaged across folds (except for computing ER) Regions ranked further (i.e. rank >37) have a perfectly null contribution to the model. ER stands for Expected Ranking. The region size is displayed in voxels. The weights were averaged across folds (except for computing ER). A region label ending in 'L' (resp. 'R') means left (resp. right) hemisphere region Table 14 displays the regions with non-null contributions to the MKL-Brodmann model, along with their contribution d m , size (in voxels) and expected ranking ER, across folds.    Regions ranked further (i.e. rank >66) have a perfectly null contribution to the model. ER stands for Expected Ranking. The region size is displayed in voxels. The weights were averaged across folds (except for computing ER)  Regions ranked further (i.e. rank >73) have a perfectly null contribution to the model. ER stands for Expected Ranking. The region size is displayed in voxels. The weights were averaged across folds (except for computing ER). A region label ending in 'L' (resp. 'R') means left (resp. right) hemisphere region  Regions ranked further (i.e. rank >46) have a perfectly null contribution to the model. ER stands for Expected Ranking. The region size is displayed in voxels. The weights were averaged across folds (except for computing ER) Regions ranked further (i.e. rank >85) have a perfectly null contribution to the model. ER stands for Expected Ranking. The region size is displayed in voxels. The weights were averaged across folds (except for computing ER)

Regions Selected in MKL-HCP
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.