Zaruhi Alaverdyan

Julien Jung

Romain Bouet

Carole Lartizien

Regularized siamese neural network for unsupervised outlier detection on brain multiparametric magnetic resonance imaging: application to epilepsy lesion screening

Keywords: Regularized Siamese network, Wasserstein autoencoder, Unsupervised representation learning, Brain lesions, Anomaly detection, Deep Learning

come     Regularized siamese neural network for unsupervised outlier detection on brain multiparametric magnetic resonance imaging: Application to epilepsy lesion screening

Introduction

A number of neurological pathologies are characterized by small lesions of various shapes, localizations and spatial patterns. This encompasses small vessel disease (SVD) [START_REF] Wardlaw | Neuroimaging standards for research into small vessel disease and its contribution to ageing and neurodegeneration[END_REF]), intracranial carotid artery calcification (ICAC) [START_REF] Bos | Intracranial carotid artery atherosclerosis and the risk of stroke in whites: the rotterdam study[END_REF][START_REF] Bos | Atherosclerotic calcification relates to cognitive function and to brain changes on magnetic resonance imaging[END_REF]) which both may lead to stroke as well as cognitive impairment or dementia, multiple sclerosis (MS) [START_REF] Filippi | Mri criteria for the diagnosis of multiple sclerosis: Magnims consensus guidelines[END_REF]) and medically refractory epilepsy, which is often associated with malformations of cortical development (MCD), such as heterotopia or focal cortical dysplasia (FCD) [START_REF] Guerrini | Epilepsy and malformations of the cerebral cortex[END_REF]). For all these pathologies, the characterization of lesion profiles is crucial to perform an early diagnostic as well as define and monitor the optimal therapeutic strategy.

Neuroimaging techniques, especially those based on multiparametric magnetic resonance imaging (MRI), have been exploited to detect, characterize and monitor this type of pathologies in a non-invasive manner. The detection of small and subtle brain lesions during standard visual examinations, however, remains challenging. As an example, recent retrospective studies involving surgical epilepsy patients indicate that 30% to 80% of lesions, depending on their type, go undetected during routine MRI exams [START_REF] Bernasconi | MRI-negative epilepsy: evaluation and surgical management[END_REF]), dramatically impacting the success of resective surgery.

Considering the impact the detection of subtle brain lesions on MR imaging has on the diagnostic and therapeutic orientation, the development of automated neuroimaging analysis tools, providing probabilistic maps of the suspicious regions, may be of valuable help to neurologists during the diagnostic phase or therapeutic monitoring. The automatic detection of subtle brain pathologies, including MS, SVD and epilepsy lesions, on multiparametric MR imaging has been increasingly addressed over the past ten years. The advances in machine learning and, more recently, deep learning have further motivated new studies in this domain. A vast majority of the existing works was performed in the supervised learning framework, thus requiring large and accurately annotated data sets. However, acquiring a data set adequately representing the heterogeneity of the pathology at hand is a major issue, more so for the lesions that may be located anywhere in the organ of interest, have various shapes, size and texture, or even can not be visually identified on the images.

In this study, we tackle the problem of subtle brain lesion detection as an outlier detection problem, to both encompass the lack of annotated pathological cases and the class imbalance problem impacting the performance of supervised learning models. We consider this challenging detection problem from the perspective of unsupervised learning combining data-driven deep feature extraction and outlier detection in a single framework. We specifically target the detection of cryptogenic epilepsy lesions in multiparametric MRI, focusing on MRI-negative lesions, meaning that these lesions were not visually identified by clinicians on the MR scans, as illustrated on figure 1.

This work builds on the method that we proposed in El Azami et al. ( 2016)

Figure 1: Examples of three patients with confirmed epileptogenic lesions on FLAIR (left) and T1-weighted MRI scans, respectively. The ground truth is annotated in red circles. The first two lesions were detected during a visual analysis (MRI-positive); the third patient has an anomaly in the hippocampus that was not detected on MRI (MRI-negative).

to detect FCD and heterotopia lesions on T1-weighted (T1w) MR images. This pipeline combines a hand-crafted feature extraction mechanism targeting typical FCD lesions with an outlier detection model based on one-class support vector machine (oc-SVM). We generalize this framework to detect a wider range of subtle brain anomalies by exploiting unsupervised deep learning architectures as feature extraction mechanisms.

Our contributions in this work are summarized below:

1. We define a novel unsupervised representation model of normal brain patterns based on a siamese network composed of stacked convolutional autoencoders as subnetworks and a loss function tailored to the context of outlier detection.

2. This unsupervised deep representation is coupled to a one class SVM model to derive a brain anomaly detection model highlighting suspicious regions in 3D brain MRI at voxel level. 3. This pipeline is applied to the detection of MRI negative lesions in multiparametric MRI (T1w and FLAIR) with competitive performances with regards to state-of-the art methods.

A proof of concept of our anomaly detection pipeline was presented in Alaverdyan et al. (2018a,b). In this work, we build on this preliminary study, first, by providing a detailed description of the model and performing an extensive performance analysis based on a clinical dataset including 75 healthy subjects and 21 epilepsy cases (including 18 MRI negative patients). We then compare our unsupervised deep latent representation model based on the early fusion of T1w and FLAIR MRI with other unsupervised deep models including Wasserstein autoencoders as well as with intermediate fusion strategies based on multiple kernel learning. Finally, we put our results in perspective with those achieved by the most recent state of the art methods.

The paper is organized as follows. We start by reviewing the existing research in related domains before describing the components of our automated pipeline in Section 3. This includes the presentation of our unsupervised representation learning model based on a regularized siamese network. Section 4 describes the different experiments conducted to evaluate the performance of our brain anomaly detection models and compare it with alternate deep unsupervised architectures. Finally, Section 5 illustrates the obtained results, followed by a discussion in Section 6.

Related work

2.1. Automated subtle brain lesion detection with deep supervised or weakly supervised architectures The recent advances in deep learning and promising performances achieved in the domain of medical image segmentation [START_REF] Litjens | A survey on deep learning in medical image analysis[END_REF]; [START_REF] Kamnitsas | Efficient multi-scale 3d cnn with fully connected crf for accurate brain lesion segmentation[END_REF]) motivated new studies, tackling the detection of subtle brain pathologies as a supervised segmentation task, e.g. for MS lesion segmentation [START_REF] Hashemi | Asymmetric loss functions and deep denselyconnected networks for highly-imbalanced medical image segmentation: Application to multiple sclerosis lesion detection[END_REF][START_REF] Valverde | Improving automated multiple sclerosis lesion segmentation with a cascaded 3d convolutional neural network approach[END_REF][START_REF] Brosch | Deep 3d convolutional encoder networks with shortcuts for multiscale feature integration applied to multiple sclerosis lesion segmentation[END_REF], [START_REF] Havaei | Hemis: Heteromodal image segmentation[END_REF]), for ICAC [START_REF] Bortsova | Segmentation of intracranial arterial calcification with deeply supervised residual dropout networks[END_REF]). Some recent studies introduced specific losses accounting for the class imbalance [START_REF] Hashemi | Asymmetric loss functions and deep denselyconnected networks for highly-imbalanced medical image segmentation: Application to multiple sclerosis lesion detection[END_REF]; [START_REF] Sudre | Generalised dice overlap as a deep learning loss function for highly unbalanced segmentations[END_REF]). Although these deep architectures achieve impressive results compared to the state-of-the art performance, even for the segmentation of small lesions, they still require large voxel-wise annotated data sets for training. Moreover, as far as we know, such supervised architectures are tailored to segmentation of lesions that can be visually detected on the image; their performance on barely visible lesions has not been reported yet. Some authors recently proposed to formulate segmentation tasks in semi-or weakly-supervised deep settings which allows to reduce or even bypass the need of voxel-level annotated data sets [START_REF] Cheplygina | Not-so-supervised: a survey of semi-supervised, multi-instance, and transfer learning in medical image analysis[END_REF]). A few attempts have been made to apply these methods to the problem of subtle lesion detection. [START_REF] Baur | Semi-supervised deep learning for fully convolutional networks[END_REF] introduced a framework for MS lesion segmentation in multi-parametric MRI coupling a standard U-Net architecture [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]), trained on labeled data, with manifold embedding accounting for unlabeled data. This model shows promising performance, except for the detection of very small lesions. In [START_REF] Dubost | Gp-unet: Lesion detection from weak labels with a 3d regression network[END_REF], the authors exploit weak labels (the number of lesions in a scan) in a U-Net like architecture to segment SVD lesions in the basal ganglia based on PD-weighted MRI scans. This framework allows achieving a sensitivity 20% higher than the state-of-the art methods. The method, however, requires to input the number of lesions, thus assuming that the detection task can be performed easily. Such methods are indeed very promising to perform fastidious segmentation tasks with weak supervision. The existing methods, however, have not been designed or optimally tuned to perform challenging detection tasks.

Deep unsupervised anomaly detection problem

Another recent tendency specifically casts the lesion detection problem as an anomaly detection task. Anomaly detection, also referred to as outlier de-tection, consists in identifying the observations that seem to be drawn from a different distribution than the one generating the normal examples. Over the recent years, the challenging topic of outlier detection has been studied extensively in different application domains [START_REF] Chandola | Anomaly detection: A survey[END_REF] and has also been recently addressed from the perspective of deep learning [START_REF] Kiran | An overview of deep learning based methods for unsupervised and semi-supervised anomaly detection in videos[END_REF]).

One category of deep anomaly detection methods is based on projecting the data samples to a low dimensional manifold and then mapping them back to the original space through a reconstruction of the original data points. The reconstruction error is later used to distinguish anomalies. In different computer vision studies, the reconstruction was obtained via various deep architectures such as autoencoders (AE) in Munawar et al. (2017a), variational autoencoder (VAE) in [START_REF] An | Variational autoencoder based anomaly detection using reconstruction probability[END_REF], long short memory networks (LSTM) in Munawar et al. (2017b) or generative adversarial networks (GAN) in [START_REF] Zenati | Efficient gan-based anomaly detection[END_REF]. In the medical imaging domain, [START_REF] Schlegl | Unsupervised anomaly detection with generative adversarial networks to guide marker discovery[END_REF] were the first to propose a GAN-based architecture, trained on normal samples only, associated to an anomalous score function composed of the GAN reconstruction and discrimination losses. This architecture was applied to the detection of fluid in high resolution clinical optical coherence tomography of the retina. Other architectures such as autoencoders [START_REF] Pawlowski | Unsupervised lesion detection in brain ct using bayesian convolutional autoencoders[END_REF]) or adversarial autoencoders (AAE) [START_REF] Chen | Unsupervised detection of lesions in brain mri using constrained adversarial auto-encoders[END_REF]) were adapted to the brain tumor segmentation problem.

Another category of anomaly detection methods seeks to learn a discriminative boundary around the normal training instances [START_REF] Chandola | Anomaly detection: A survey[END_REF]). Any test instance that does not fall within the learnt boundary is declared anomalous. From the deep learning perspective, this amounts to first learning latent representations of normal samples with a deep unsupervised network, similar tot the first category of anomaly detection methods cited above, and then feeding the learned representations to a one-class classification algorithm in order to estimate the boundaries of the normal examples, as proposed by [START_REF] Erfani | Highdimensional and large-scale anomaly detection using a linear one-class svm with deep learning[END_REF]. As far as we know, we are the first to build on this approach in the medical image analysis domain and for the challenging task of subtle epilepsy detection (Alaverdyan et al. (2018a,b)).

State-of-the-art automated detection systems for epilepsy lesion detection

Medically refractory epilepsy is often associated with malformations of cortical development (MCD) encompassing a large diversity of lesion types [START_REF] Barkovich | A developmental and genetic classification for malformations of cortical development: update 2012[END_REF]). MCD include focal cortical dysplasia (FCD) categorized as histological subtypes I, II and III as well as focal, band-shaped or subependymal subcortical heterotopia characterized by the presence of neurons located deep in the white matter. These lesions may also be associated to gliosis and hippocampus sclerosis. Recent retrospective studies involving surgical epilepsy patients indicate that up to 33% with typical FCD type II lesions and 87% with FCD type I lesions go undetected during routine MRI exams [START_REF] Bernasconi | MRI-negative epilepsy: evaluation and surgical management[END_REF]). Similarly, subtle heterotopia may only become apparent after MRI post-processing (Huppertz et al. (2005)). The success rate of surgery is around 70% [START_REF] Wiebe | A randomized, controlled trial of surgery for temporal-lobe epilepsy[END_REF]; [START_REF] Keller | Persistent seizures following left temporal lobe surgery are associated with posterior and bilateral structural and functional brain abnormalities[END_REF]; [START_REF] Bien | Trends in presurgical evaluation and surgical treatment of epilepsy at one centre from 1988-2009[END_REF]) but is significantly lower in MRI-negative patients than for MRI-positive patients [START_REF] Alarcon | Is it worth pursuing surgery for epilepsy in patients with normal neuroimaging[END_REF]; [START_REF] Bell | Epilepsy surgery outcomes in temporal lobe epilepsy with a normal mri[END_REF]; [START_REF] Bien | Characteristics and surgical outcomes of patients with refractory magnetic resonance imaging-negative epilepsies[END_REF]).

Over the recent years, various automated detection systems have been proposed for the challenging task of epilepsy lesion detection [START_REF] Kini | Computational analysis in epilepsy neuroimaging: A survey of features and methods[END_REF]). The vast majority of those systems specifically target FCD type II lesions and address this problem from a supervised machine learning perspective combining standard classification algorithms (e;g SVM, decision trees) with manually engineered features that have been shown to correlate with the appearance of these lesions on MR scans, including features derived from surface based morphometry (SBM) (Hong et al. (2014a); [START_REF] Ahmed | Cortical feature analysis and machine learning improves detection of "mrinegative" focal cortical dysplasia[END_REF]; [START_REF] Gill | Automated detection of epileptogenic cortical malformations using multimodal mri[END_REF]). More recently, [START_REF] Gill | Deep convolutional networks for automated detection of epileptogenic brain malformations[END_REF] proposed to harness deep learning by designing an architecture composed of two convolutional neural networks trained on T1w and FLAIR patches extracted from the gray matter area in the brain. They demonstrated an improved detection performance compared to the method proposed in [START_REF] Gill | Automated detection of epileptogenic cortical malformations using multimodal mri[END_REF] combining a number of SBM, intensity and gradient features with an ensemble of RUSBoosted decision trees. All these approaches require a delineation of lesional zones which is not accurate in MRI-negative patients. For these patients, the ground truth comes from the post-surgical scans containing the resected zone which also includes non-lesional voxels. In [START_REF] Ahmed | Cortical feature analysis and machine learning improves detection of "mrinegative" focal cortical dysplasia[END_REF], the authors showed that manually reducing the resection masks for MRI-negative patients to correct the label noise resulted in a detection rate of 58% while more "generous" annotations achieved only 12%.

To bypass the difficulty of obtaining ground truth annotations, other approaches propose to compare patients to a cohort of normal subjects in order to discriminate the lesions [START_REF] Thesen | Detection of epileptogenic cortical malformations with surface-based mri morphometry[END_REF][START_REF] Srivastava | Feature-based statistical analysis of structural mr data for automatic detection of focal cortical dysplastic lesions[END_REF][START_REF] Bruggemann | Voxel-based morphometry in the detection of dysplasia and neoplasia in childhood epilepsy: combined grey/white matter analysis augments detection[END_REF]). Voxel-based morphometry (VBM) analysis builds on this approach by computing mass univariate general linear models (GLM) at the voxel level based on a feature map encoding typical FCD lesion patterns such as cortical thickening and blurring of the GM/WM interface [START_REF] Srivastava | Feature-based statistical analysis of structural mr data for automatic detection of focal cortical dysplastic lesions[END_REF][START_REF] Bruggemann | Voxel-based morphometry in the detection of dysplasia and neoplasia in childhood epilepsy: combined grey/white matter analysis augments detection[END_REF]. In another unsupervised setting, [START_REF] Ahmed | Decrypting "cryptogenic" epilepsy: Semi-supervised hierarchical conditional random fields for detecting cortical lesions in mri-negative patients[END_REF] evaluated the performance of hierarchical conditional random fields on one among four SBM features or their combination, focusing on cryptogenic epilepsy. In El Azami et al. ( 2016), we extracted textural features (Huppertz et al. (2005); Wagner et al. (2011)) to learn a one class SVM (oc-SVM) model per voxel and later identify the lesions as the clusters of voxels with the most negative scores output by oc-SVMs. The reasoning behind the choice of building location specific (e.g. voxel based) classification models is that it should enable detecting subtle pattern variations induced by the presence of the lesion, unlike global models pooling all voxels.

Limitations of the current computational models for epilepsy detection have been recently pointed out in survey by [START_REF] Kini | Computational analysis in epilepsy neuroimaging: A survey of features and methods[END_REF]. In particular, the au-thors emphasize the need to acknowledge all types of epilepsy lesions, and recommand to circumvent the lack of massive databases representing the variability of all pathological cases by considering the development of neurocomputational models, first, from the perspective of unsupervised outlier detection and, second, by combining complementary information provided by multimodality imaging to gain in specificity.

Method

In this work, we build on the anomaly detection frameworks proposed by [START_REF] Erfani | Highdimensional and large-scale anomaly detection using a linear one-class svm with deep learning[END_REF] and El Azami et al. ( 2016) by combining unsupervised feature learning and robust one class classification with oc-SVM for anomaly detection. In the following, we present the general pipeline of our anomaly detection system, then detail the different elements of this pipeline, especially emphasizing the novel architecture of unsupervised representation model based on siamese networks.

General pipeline

In the proposed general anomaly detection system, each voxel is characterized by a latent representation vector extracted from an unsupervised deep network. For each voxel in the brain, its normality is modeled with a oc-SVM classifier in the latent representation space. Eventually, for an unseen patient, abnormalities can be found as local neighborhoods of voxels found anomalous by the corresponding oc-SVM models. The general pipeline is illustrated on figure 2. It consists of two major steps • In the first step, we extract image patches of all the available volumes of the healthy control MRI dataset and learn a latent representation of each patch with a deep unsupervised network. Once this step is performed, the central voxel of each patch extracted from a brain volume will be associated to a latent representation yielded by this deep network.

• In the second step, we build one oc-SVM model per voxel in the latent representation space learned at the previous step. Each voxel is associated with a oc-SVM classifier, hence the number of classifiers is equal to the number of voxels in the volume of interest. For a given voxel v i , the associated oc-SVM classifier C i is trained on the matrix composed of the representations of the patches of all the subjects from the healthy control dataset centered at v i .

For a new patient, each voxel v i is first assigned a latent representation based on the deep unsupervised model; this latent representation is then matched against the corresponding classifier C i and is assigned the signed score output by the classifier. This yields a distance map D p for the given patient.

Note that the pipeline represented on fig. 2 takes one imaging modality, 3D MRI T1-W, at input, for clarity purpose. Our model actually combines T1-w an FLAIR images as input channels of the deep unsupervised architecture. Each step of the system will be explained in details in the following sections.

Regularized siamese neural network for representation learning

Rationale

Our objective is to map the original patches to a representation space that captures the global semantic information from the image whatever the brain region and where the patches belonging to different healthy subjects and centered at the same spatial localization are close with respect to a chosen metric.

Autoencoders [START_REF] Hinton | Autoencoders, minimum description length and helmholtz free energy[END_REF]) and their variations have been shown efficient for feature extraction in various contexts and applications, while siamese networks [START_REF] Bromley | Signature verification using a "siamese" time delay neural network[END_REF]; [START_REF] Chopra | Learning a similarity metric discriminatively, with application to face verification[END_REF]) perform well in learning representation space where the distance between similar and dissimilar pairs of instances can be controlled. In this study, we leverage both autoencoders and siamese network architectures in a unified framework to match our objective.

We hypothesize that such a siamese network composed of identical autoencoders as subnetworks should better capture the shared fine patterns of each patch by encoding the similarity constraint in the loss function and thus reinforce the compacity of the patch distribution in the latent space with regards to standard autoencoders. The latent representations of healthy patches centered at the same spatial location will thus be driven closer while that of an anomalous patch centered at the same location in the brain, which has never been processed by the network, will lay further in the representation space.

We also propose a novel architecture that bypasses dissimilar pairs unlike standard siamese architecture and is thus trained on similar pairs only. This choice is motivated by the fact that the notion of dissimilar pairs can not be objectively defined in our context. The inclusion of dissimilar pairs, however, is not an absolute prerequisite of the good performance of siamese architecture. [START_REF] Zheng | Triangular similarity metric learning for face verification[END_REF], indeed, showed that siamese networks learn better latent representations when trained on similar pairs only for an application to face verification. The intuition behind this result is that learning dissimilarities is challenging and not well conditioned as opposed to learning similarities. Thus, depending on the data, the two terms of the objective function related to dissimilarity and similarity may be partly contradictory and inhibit the learning of similarities during the training process.

Network Architecture and loss function

The proposed architecture is illustrated on figure 3. Our regularized siamese neural network (rSN) consists of two identical (same architecture, shared parameters) subnetworks -stacked convolutional autoencoders (CAE) with K hidden layers and a cost module. The input x of a CAE is first encoded to a middlelayer representation by a series of convolutional and max-pooling operations and later decoded with a series of deconvolutions and up-poolings to produce a reconstruction x of the input. A convolutional layer l is composed of N l kernels and biases and can be expressed as

H m l = f (W m l-1 * H l-1 + b m l-1 )
where H m l is the m-th feature map of layer l, W m l is the kernel matrix associated with H m l and b m l is its bias, f is an activation function (usually non-linear). * denotes the convolution operation. The siamese network receives a pair of patches (x 1 , x 2 ) at input, then each patch is propagated through the corresponding subnetwork yielding representations z t ∈ Z, t = (1, 2) in the middle layer which are then passed to the loss function L rSN below. Our loss function is designed to maximize the cosine similarity between z 1 and z 2 . To compensate for the lack of dissimilar pairs, we propose to add a regularizing term consisting of the mean squared error between the input patches and their reconstructions output by the subnetworks. Without a proper regularization term, the loss function could be driven to 0 by mapping all the patches to a constant value. The proposed loss function for a single pair hence is:

L rSN (x 1 , x 2 ; Θ) = 2 t=1 ||x t -xt || 2 2 -α • cos(z 1 , z 2 ) (1)
where xt is the reconstructed output of subnetwork t of the patch x t while z t is its (vectorized) representation in the middle layer and α is an hyperparameter that controls the tradeoff between the two terms. Θ represents the parameter set.

Voxel-level outlier detection with oc-SVM

A one class support vector machine (oc-SVM) [START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF]) is an outlier detection method based on the SVM algorithm assigning labels y i ∈ {-1, 1} to two distinct classes of objects, based on n training samples (z i , y i ) ∈ X from the negative class only. The training examples are first mapped to a higher dimensional space via a feature map φ associated with a kernel K such that K(z i , z j ) = φ(z i ), φ(z j ) . The corresponding optimization problem is the following:

min w,ρ,ξi 1 2 ||w|| 2 -ρ + 1 νn n i=1 ξ i subject to w • φ(z i ) ≥ ρ -ξ i , ξ i ≥ 0, i ∈ [1, n] (2)
where ξ i -s are slack variables relaxing the inequality constraints so as to account for the non-separable classes, w and ρ define the separating hyperplane, ν is a parameter that sets a boundary to the fraction of allowed outliers. The decision function, for an example z, is w • φ(z) -ρ. This decision function contributes to the signed score output by a oc-SVM model. In a typical scenario, examples with negatives scores would be considered outliers.

In the scenario depicted on figure 2, a given voxel v i is associated with a oc-SVM classifier C i trained on the matrix M i = [z i1 , ..., z in ] where z ij is the feature vector corresponding to the patch centered at v i of subject j and n is the number of subjects. During the test phase, each voxel v i of a given patient p is matched against the corresponding classifier C i and is assigned the signed score output by the classifier. This yields a distance map D p for the given patient.

Post-processing

For a given patient p, the distance map D p outputted at the previous step is then post-processed to obtain the final anomaly detection map. A 3-step post-processing is proposed as follows.

The first step consists in normalizing the distance maps with respect to the intra-subject spatial variability. For that purpose, the distance maps of the control subjects are computed by performing a k-fold evaluation of the training set (i.e. for each fold of normal subjects, the distance maps are obtained with oc-SVMs trained on the remaining subjects). These maps are used to estimate the standard deviation of the normal subjects' distance distribution at voxellevel. For a given patient p, a new map Dp is computed by a voxel-wise division of the output distance map D p over the estimated standard deviations. The final distance map F p is then derived by averaging D p and Dp i.e.

F p = 1 2 ( D p max(abs(D p )) + Dp max(abs( Dp ))
).

(3)

The reason behind the additional term is that some zones in the brain have more intra-subject variability than others and therefore are more likely to be considered as anomalies. By weighing them by the standard deviation, the score maps account for this effect.

The second step consists in thresholding the F p map to produce a cluster map. To this end, all the voxel score values of F p are pooled together into a histogram which is then approximated by a non-parametric distribution using a kernel density estimator [START_REF] Bowman | Applied smoothing techniques for data analysis: the kernel approach with S-Plus illustrations[END_REF]). The approximated patient distance score distribution is then thresholded at some pre-chosen value and a 26-connectivity rule is applied to identify the connected components. These components are referred to as clusters. By varying the threshold, the number of clusters can be controlled according to a clinician's needs. An example of F p output score map, as well as the cluster maps obtained with different thresholds are shown on fig. 4. Clusters smaller than a certain size may also be discarded.

The third step consists in ranking the detected clusters. We use the following ranking criterion to assign a rank to a cluster c i , inspired from Ahmed et al. ( 2016)

rank(c i ) ∼ ω * score(c i ) min j score(c j ) + (1 -ω) * size(c i ) max j size(c j ) ( 4
)
where score(c i ) is the average of the voxel scores in the cluster and size(c i ) is the number of voxels in the cluster and ω is a parameter weighing the relative contribution of the cluster size and score. Such a ranking favors large clusters with the most negative average score. Using this ranking, we keep the top n detections and discard the rest.

Experimental settings

Dataset description and pre-processing

The study was approved by our institutional review board with approval numbers 2012-A00516-37 and 2014-019 B and a written consent was obtained for all participants.

Study group

The data set considered in this study consists in a training set of healthy subjects and a test set of epilepsy patients. Patient group: The test group consists of 21 patients who had been admitted to the Neurological Hospital of Lyon and diagnosed with medically intractable epilepsy. The age of the patients varies between 17 and 47 years, with a median of 29. As a part of the pre-surgical evaluation, they all had T1-weighted (T1w) and FLAIR MRI sequences as detailed below. Additionally, the patients underwent intracranial EEG exams in order to localize the origin of seizures. Healthy control group: The training data set consists of 75 healthy individuals aged between 20 and 66 years. All the subjects had similar T1w and FLAIR MRI sequences as the patient group. Sixty-nine of these control subjects were used as the training datasets while the remaining 6 subjects were used as a validation dataset.

MRI acquisition and pre-processing

All the healthy controls and patients had 3D anatomical T1w brain MRI sequences (TR/TE 2400/3.55; 160 sagittal slices of 192 x 192 1.2mm cubic voxels) and FLAIR MRI sequences (176 slices of 196 x 256 1.2mm cubic voxels) on a 1.5 T Sonata scanner (Siemens Healthcare, Erlangen, Germany). The pre-processing was done using the SPM8 software [START_REF] Ashburner | Computational anatomy with the spm software[END_REF]). The T1-weighted MRI volumes were first processed with the unified segmentation algorithm (UniSeg) [START_REF] Ashburner | Unified segmentation[END_REF]) implemented in SPM (using the default parameter values) that performs tissue segmentation (white/grey matter, cerebrospinal fluid), correction for magnetic field inhomogeneities and spatial normalization. All the 3D MR volumes were normalized to the standard brain template of the Montreal Neurological Institute (MNI) [START_REF] Mazziotta | A probabilistic atlas and reference system for the human brain: International consortium for brain mapping (icbm)[END_REF]) with a voxel size of 1 x 1 x 1 mm. Further, FLAIR images were rigidly co-registered to the individual T1w MR images. Next, the transformation from the subjects' native space to the MNI space, produced by the UniSeg algorithm, was applied on the co-registered FLAIR images in order to normalize them to the MNI space as well.

We excluded the brain regions (the cerebellum and brain stem) that are not susceptible to epilepsy using a masking image in the MNI space derived from the Hammersmith maximum probability atlas described in [START_REF] Hammers | Three-dimensional maximum probability atlas of the human brain, with particular reference to the temporal lobe[END_REF].

After the elimination of the corresponding voxels, the number of remaining voxels adds up to around 1.5 million per volume. We removed top 1% intensities and scaled the images between 0 and 1 at image level before feeding the patches to the deep architectures.

Patient lesion annotations

The patient group was selected in collaboration with the expert neurologists of the Lyon Neurological University Hospital to sample highly difficult detection cases so as to challenge the performance of the proposed automated detection system. Most of them had a resective surgery and became seizure free at most 6 months after the surgery. A few patients had successful thermocoagulation instead (patients A -, G + ). The positive outcome of resective surgery or thermocoagulation is considered as the ground truth to establish the epileptogenic lesion localization. The ground truth annotations used in the performance evaluation were thus obtained by outlining the visible zones of the MRI-positive patients (D + , G + and R + ), and by combining the information from the intracranial EEG, post-surgical or post-thermocoagulation MR images and the resected zones for MRI-negative patients. Note that the MRI-negative patients considered in this study remained MRI-negative after a retrospective inspection of the T1w and FLAIR scans. For patients undergoing surgery, the histopathological analysis of the resected tissue indicated FCD type II in 3 patients (D + , G + and Q -) and FCD type III in 3 patients (C -, P -and S -). The histopathological analysis was inconclusive for the remaining patients, which is in accordance with the statistics reported in [START_REF] Bernasconi | Advances in mri for "cryptogenic" epilepsies[END_REF] and underlines our objective to detect lesions with unknown signatures. A full description of the lesion types and localizations is provided in Table 1.

Implementation details of the brain anomaly detection model

The general architecture of the brain anomaly detection model depicted on figure 2 was adapted to the specific application of epilepsy lesion detection in multiparametric T1-w and FLAIR MRI. Hereafter, we provide the detailed description of the different elements of the pipeline. The entire pipeline was implemented in Python, with Theano/Keras libraries for the representation learning stage and a LibSVM wrapper in scikit-learn [START_REF] Pedregosa | Scikitlearn: Machine learning in Python[END_REF]) for one class SVM.

Unsupervised representation learning with rSN.

The proposed rSN consists of two identical stacked convolutional autoencoders as depicted on figure 3 and on the top left part of figure 6. The specific architectures of the encoder (E ) and decoder (or generator) (G) elements are shown on figure 5. The input of the encoder consists of the patches of each modality combined as channels. For each couple of patches of a given subject, a random 'similar pair' is selected among the identically located patches of other subjects yielding in total around 3.5 million pairs of patches for training. The input patch size was set to 15x15 after a number of tested configurations and is justified by the subtle nature and size of epilepsy lesions. Indeed, larger patch sizes (e.g. 30x30) were not successful at detecting subtle lesions. The encoder and decoder are composed of two convolutional layers with kernel size of 3x3 and a stride of 1. A max-pooling operation is performed only at the first and last layer of the encoding and decoding part respectively. As shown on figure 5, the middle layer has 16 feature maps of 2x2 which, when flattened, yields a 64-dimensional vector. This dimension compromises the need to capture subtle patterns at the patch level while conforming to the size of the training database (69 subjects), to avoid over-fitting of the oc-SVM models. We used ReLU activation function in all the layers except the last one where the sigmoid function is used. The α parameter in loss (1) was set to 0 during the first 10 epochs, then grew linearly for 15 epochs until it reached some α max value and then plateaued for 5 more epochs. The Adam optimizer was used with the learning rate set to 0.001.

outlier detection with oc-SVM

We used oc-SVM classifiers with RBF kernel. The γ kernel parameter was derived for each voxel v i individually by setting it to the median of the standardized euclidean pairwise distances of the corresponding matrix M i (see section 3.3) as in [START_REF] Caputo | Appearance-based object recognition using svms: which kernel should i use?[END_REF]. Varying the parameter ν corresponding to the upper bound on the fraction of permitted outliers (see eq. 2) did not significantly impact the results; the fraction of the outliers is indeed controlled in the post-processing step with the threshold applied on the distance map F p . Thus, ν was set to 0.03 for all voxels and all scenarios.

Post-processing

In the first post-processing step, we performed a 10-fold cross validation to compute the output maps of the normal population that served to derive the output map D in eq. 3. In the second step of the post-processing, we empirically set the threshold of the output score maps to the value that resulted in at most 10 clusters. We indeed observed among the output maps of the normal population that larger values typically produce a small number of very large clusters. The minimal cluster size was set to 82 voxels corresponding to the expected cluster size calculated with the SPM analysis of the T1w MRI data (see supplementary file). The size of the majority of the detected clusters varies between 500 and 3000 voxels, this threshold therefore does not affect the performance in any significant way. In the third post-processing step, the parameter ω weighing the relative contribution of the cluster size and score was set to 0.5 to equally account for the cluster average score and size.

Comparison with intermediate fusion strategy of T1-w and FLAIR modalities

We explore an alternative strategy to the early fusion of T1w and FLAIR MRI modalities as input channels of the siamese network. This consists in performing intermediate fusion by training individual networks for each modality and combining the learned representations with a multiple kernel algorithm. We propose to use the slimSimpleMKL algorithm [START_REF] Loosli | Using svdd in simplemkl for 3d-shapes filtering[END_REF]) to leverage the original formulation of the simpleMKL algorithm [START_REF] Rakotomamonjy | Simplemkl[END_REF]) that extends the oc-SVM formalism to fit the multiple kernel paradigm [START_REF] Bach | Multiple kernel learning, conic duality, and the smo algorithm[END_REF]; [START_REF] Sonnenburg | A general and efficient multiple kernel learning algorithm[END_REF]. By controlling the number of support vectors with a tradeoff parameter λ, tight normality bounds can be achieved with SlimSimpleMKL which, in turn, can lead to an improved performance of the outlier detection model. Two separate architectures referred to as T1 rSN and FLAIR rSN were thus trained with the T1-w and FLAIR healthy subject dataset, respectively. Both architectures used similar encoder and decoder architectures as those depicted on figure 5 except that the number of input channel was set to 1. The intermediate fusion with slim-SimpleMKL was performed in Matlab, using the implementation provided by [START_REF] Loosli | Using svdd in simplemkl for 3d-shapes filtering[END_REF].

Comparison with other deep unsupervised representation models

We compare the proposed siamese model with two alternative unsupervised representation models, stacked convolutional autoencoders (CAE) and Wasserstein autoencoders (WAE) [START_REF] Tolstikhin | Wasserstein Auto-Encoders[END_REF]). Figure 6 illustrates the general setup. In our comparison, the encoder E and decoder G of CAE and WAE have the same architecture as the ones of the rSNN subnetworks shown on figure 5.

The stacked convolutional autoencoder (CAE) serves as a baseline performance model to evaluate the potential added value of the proposed rSN model and its ability to capture a finer representation of the normal brain based on the proposed loss L rSN .

Wasserstein autoencoders have been recently introduced as generative models combining the best properties of Wasserstein GANs and Variational autoencoders [START_REF] Tolstikhin | Wasserstein Auto-Encoders[END_REF]). They consist of three components: an encoder E, a decoder G and an adversary network D that tries to distinguish the prior distribution of the latent code P Z from the latent distribution Q Z produced by the encoder. The resulting loss function can be expressed as

L W AE (X; Θ W AE ) = 1 N N i=1 c(x i , xi ) + β • D Z (P z , Q z ) (5)
where D Z measures the discrepancy between a given distribution P z and Q z for the dataset X = {x i } 1,..,N and c measures the reconstruction error. β is a coefficient that controls the tradeoff between the two terms and Θ W AE denotes the parameter set. The generic form of the WAE loss allows different reconstruction error functions and regularizers. We used the standard reconstruction error

c(x i , xi ) = ||x -xi || 2
2 and the Jenssen-Shanon divergence as D Z , estimated with a discriminator. As for the proposed rSN model, we hypothesize that this architecture should capture a finer representation than the standard CAE by enforcing the learned representations to follow the prior distribution. In our comparison, the WAE discriminator D consists of four fully connected layers of dimension 128, 128, 64 and 1 with LeakyReLU as activation (with scale 0.02 for negative input values). We varied the parameter β in the L W AE loss (eq. 5) among the following values -1,5,10,20 and 100.

Performance evaluation

We chose typical metrics applied in lesion detection tasks to evaluate the performance of our brain anomaly detection model. We define the number of true positive (TP) and false positive (FP) detections at cluster-level i.e. we consider that a cluster is true positive when it overlaps with the ground truth annotation and false positive (FP) otherwise. As such, we quantify the sensitivity of the system as the percentage of the patients whose lesions were correctly detected. We couple the system sensitivity with the average number of FPs per patient to plot a fROC curve [START_REF] Bunch | A freeresponse approach to the measurement and characterization of radiographicobserver performance[END_REF]). We also perform a qualitative visual analysis of the detected cluster maps and provide a quantitative evaluation at the patient level.

Results

We first evaluate the performance of the proposed brain anomaly detection model depicted on figure 2. We then compare its performance to those achieved with variants of the same global architecture. The brain anomaly detection model is trained on 69 healthy subjects and is evaluated on 21 patients with confirmed epilepsy lesions described in section 4.1. In this section, we report the performance of our brain anomaly detection model based on deep unsupervised representation extracted from a multichannel T1/FLAIR regularized siamese network and associated to a oc-SVM classification and ad-hoc post-processing.

Performance of our brain anomaly detection model

Figure 8a reports the fROC curves computed on the series of 21 patients considering two values of α1 in loss L rSN (see eq (1)) -0.25 and 0.5. Table 1 reports performance at patient level for the best configuration corresponding to α = 0.5.

The system achieves an overall sensitivity of 62% (13/21) with a mean rank of the TP detections of 3.6 on the 21 patients. The corresponding performance on the 18 M RI -patients is 61% (11/18) with a mean TP rank of 4.1. The system fails to detect the lesions of unknown type in 5 patients and 2 patients with hippocampal sclerosis.

Figure 7 visualizes example cluster maps output by the brain anomaly detection model, while table 1 reports its performance at patient level. All the detected clusters, with a rank higher or equal to that of the true positive detection, are shown on the figure. The highest ranked cluster, corresponding to the most suspicious detection, is shown in light yellow while the lowest ranked clusters are depicted in darker colors ranging from light orange to dark brown. The symbol in table 1 indicates that the lesion was correctly detected by the system while the rank of the true detection is given inside parentheses. The first column on figure 7 corresponds to the MRI-positive patient D + where the highest ranked detected cluster matches the lesion localization. Table 1 confirms that the two MRI-positive patients (D + and G + ) with confirmed FCD type II lesions were detected with the highest rank. Our system allows detecting very small and subtle lesions with unknown signatures in patients N -and B -. These lesions are detected with the highest rank as indicated in table 1. Patient E -illustrates a very difficult case of a small lesion of an unknown type. In such difficult cases, the anomalies accompanying the epilepsy lesions are so subtle that they appear together with other anatomical outliers. For this patient, the lesion was ranked 8th.

For comparison, we also report the performance of the two anomaly detection models based on rSN representations learned separately on T1w and FLAIR MR images, respectively. These two models are referred to as T1 rSN and FLAIR rSN models in the following. Fig. 8a shows the corresponding fROC curves for two values of the α coefficient. Table 1 summarizes their performance at patient level. Results reported on figure 8a and table 1 clearly demonstrate that our model outperforms those based on features learned from one imaging modality. The best configuration of T1 rSN model (corresponding to α = 0.25) allowed to achieve 39% sensitivity (7/18) with a mean TP rank of 3.2 ± 2.5, for MRI-negative lesions. The best performance of the FLAIR rSN model was also reached for α = 0.25, leading to 50% sensitivity (9/18) on MRI-negative patients with a mean TP rank of 5.0 ± 3.2. These results speak of the complementary nature of the two modalities. The deep architectures are flexible with respect to accommodating multiple modalities; hence, other modalities could be easily integrated into the architecture.

Comparison with intermediate fusion strategy of T1-w and FLAIR modalities

Fig. 8b shows comparative performances of our brain anomaly detection model with a variant of this architecture consisting in applying multiple kernel learning with slimSimpleMKL on the features learnt with independent monomodal (T1w or FLAIR) rSN models. The parameter λ controlling the number of support vectors in slimSimpleMKL was varied among the values 0, 0.05, 0.1, 0.25, and 0.5. The intermediate fusion with λ = 0.1, 0.25, 0.5 results in identical performance and, therefore, the corresponding fROC curves are seen as a single one. From this comparison, it is apparent that even a slight regularization on the number of support vectors with λ = 0.05 offers a significant improvement over the original SimpleMKL formulation (λ = 0). The sensitivity jumps from 32% to 52% for the same false positive rate (9 FPs). As it can be seen on figure 8b, the early fusion with our multichannel T1/FLAIR rSN model outperforms significantly the intermediate fusion strategy with sim-pleSimpleMKL. Combining modalities in the network training stage is thus likely Table 1: Performance of our brain anomaly detection model (T1/FLAIR rSN) and comparison with a similar architecture trained on a single imaging modality (T1 rSN or FLAIR rSN). denotes a true detection followed by its rank inside parentheses. denotes no true positive detection meaning that the lesion was not detected among the 10 highest ranked clusters detected by the model.

to better leverage the complementary information contained in T1w and FLAIR MR images while the intermediate fusion, that operates on features learned separately on each modality, may skip those properties.

Comparison with alternate unsupervised representation learning models

We next compared our brain anomaly detection model with alternate models based on the same architecture but different multichannel unsupervised represention models, namely WAE and CAE, as described in section 4.3.3. Fig. 8c reports the fROC curves for these three models. As mentioned above, the WAE parameter β was varied among 1,5,10 and 20. Again, the performance achieved with the T1/FLAIR rSN features outperforms those achieved with the features extracted with the WAE and CAE models. WAE is shown to perform better than CAE for all but one value of β. The latter confirms our hypothesis that the reconstruction error, when enhanced with a regularization term, fits better to the anomaly detection context. The WAE performance is still inferior to that of rSN which might be due to a limitation of the model itself or the experimental choice of the hyper-parameters. We can see how the performance is affected by the choice of β; the value of 20 is less successful, probably since it prioritizes too much the adversarial term.

Discussion

This study presents a novel brain anomaly detection model combining unsupervised latent representation model with one-class classification at the voxel level. We have formulated a regularized siamese network architecture that learns normal brain representations using a set of non-pathological MR volumes. The features learnt with the network do not target a specific pathology but rather allow to capture normal variability from a cohort of healthy subjects. The framework allows integrating multiple modalities and we have shown the performance gain obtained by coupling T1w and FLAIR imaging for the task of detecting subtle epilepsy lesions in MRI-negative patients. To our knowledge, this is the first attempt to extract unsupervised deep latent representation for epilepsy lesion detection. The proposed approach achieves 61% detection rate on the 18 MRI-negative patients, meaning that it detects 61% of the lesions that were not visually detected by the clinicians in these challenging cases, at the same time outperforming Wasserstein and convolutional autoencoder architectures.

As stated in Section 2.3, most current studies target the detection of FCD lesions in T1-w MRI, mainly focusing on M RI + exams and do not report the average false positive detection rate per patient but rather compute specificity as the percentage of normal controls in whom no lesion was falsely identified. Table 2 summarizes the main performance of these different models. The vast majority uses manually designed features characterizing cortical malformations based on surface based morphometry (SBM) [START_REF] Thesen | Detection of epileptogenic cortical malformations with surface-based mri morphometry[END_REF]Hong et al. (2014b); 2016) and reported its performance on our patient dataset.

As shown in table 3, a performance drop is observed when the evaluation is performed on our cohort of M RI -patients including lesions with unknown signature. Our model, on the other hand, achieves good performance for the detection of FCD Type II lesions, as reported in table 1. This comparative analysis confirms the competitive performance achieved by our brain anomaly detection model with 61% detection rate and a mean rank of 4.1 ± 2.9 on the 18 M RI -lesions. The slightly inferior performance with regards to some reported results in the literature must be considered from the perspective of the challenging evaluation conditions considered in our study, where 1) the patient cohort consists of 18 purely M RI -patients and lesions of unknown type, 2) the false positive rate is reported as an average number of FP detections estimated on the patient cohort and not on a control group.

Our brain anomaly detection model fails to identify the lesions of 8 patients. A visual analysis of the system's output for those cases seems to reveal two major reasons. For some of those patients, the raw output of the system highlighted some anomaly; however, after all the post-processing steps, those clusters were not ranked among top 10 detections. This is likely to mean that other anomalies present in the original images are considered 'anomalous' to a greater extent than the subtle epileptogenic lesions. The second category involves patients whose output score maps came out without any indication of anomaly in the zone of interest. Our future work will be aimed at analyzing more thoroughly the cases when the system fails and investigate the reasons which may lay in the approach or the input images carrying no distinct marker for the lesion at all. One potential limitation of our method is the need to register all the image volumes to perform a voxel-based analysis. El Azami et al. ( 2016) used a similar framework and showed that applying two different registration methods has no significant impact on the global performance of their oc-SVM models. In our configuration, the siamese network, receiving a pair of patches a priori centered at the same location in the brain but with possible registration inaccuracies, is likely to learn a latent representation that smooths these differences, given that most patches are registered adequately. Moreover, the raw output score maps are normalized in the post-processing step by dividing them by the estimated standard deviation among the normal population, which also penalizes the highly variable zones in the images including areas that are likely to carry registration inaccuracies. We are thus confident that our system is robust to such potential registration errors.

In this study, we chose to report the top 10 detected clusters which might seem quite a high number. The output maps, however, also provide the ranking of the detected clusters, thus allowing the radiologist to adjust the number of suspicious anomalies to visualize. As reported in the result section, the mean ranking of the true positive lesion ranges around 4, meaning that retaining the top 5 lesions may not result in a significant drop in sensitivity. As discussed previously by [START_REF] Ahmed | Decrypting "cryptogenic" epilepsy: Semi-supervised hierarchical conditional random fields for detecting cortical lesions in mri-negative patients[END_REF], we think that this ranking approach provides a natural way to focus the radiologists' attention.

There are different options to improve the diagnostic performance of the proposed system. First, some pathology-specific information could be introduced in the post-processing step, by discarding some of the detected clusters based on shape and/or localization criteria. However, automated post-processing or cascaded classification systems should be cautiously evaluated since they may result in a significant sensitivity drop as observed in [START_REF] Tan | Quantitative surface analysis of combined mri and pet enhances detection of focal cortical dysplasias[END_REF]. As shown on figure 7, some of the detected false positive clusters are indeed irregularities that can be easily discarded by a trained radiologist, especially those with a low rank. An alternative option is to move towards a semi-supervised setting by enhancing the neural network with a few 'pathological' patches that could be extracted from MRI-positive cases or after a careful analysis of retrospective MRI-negative patients, following, for instance, some ideas recently proposed in [START_REF] Shah | Abnormality detection using deep neural networks with robust quasi-norm autoencoding and semisupervised learning[END_REF]. More improvement could be achieved by accounting for the complementary information provided by different imaging modalities. T1w and FLAIR modalities, introduced as channels to our network, allowed a significant diagnostic performance gain as shown on figure 8a. We expect a further performance gain by exploiting PET imaging as recently demonstrated in [START_REF] Tan | Quantitative surface analysis of combined mri and pet enhances detection of focal cortical dysplasias[END_REF].

Finally, the proposed method is quite straightforward to implement and to apply in daily practice as the output of the system can be obtained under a couple of minutes. Moreover, the system can be applied to detect other subtle pathologies which may serve as an advantage in the clinical routine.

As encouraged by [START_REF] Kini | Computational analysis in epilepsy neuroimaging: A survey of features and methods[END_REF], we also wish to make our computational pipeline available to clinicians. Ongoing work aims at transferring the different components (image pre-processing, deep feature extraction, etc.) of our brain anomaly detection model to the Neuroimaging toolbox of the VIP platform (https://vip.creatis.insa-lyon.fr) dedicated to the simulation and processing of massive data in medical imaging [START_REF] Glatard | A virtual imaging platform for multi-modality medical image simulation[END_REF]). This web portal allows users to access the CAD system as a service and significant computing resources and storage with no required technical skills beyond the use of a web browser.

Conclusion

In this study, we presented a novel configuration of unsupervised deep architectures for anomaly detection. The proposed regularized siamese network was leveraged as a representation learning mechanism, coupled with per voxel oc-SVM models. The clinical application of the proposed framework consists of automated detection of subtle epilepsy lesions in MRI-negative patients on T1-weighted and FLAIR MRI sequences. The approach achieved a sensitivity of 61% for 9 false detections per patient on pure MRI-negative patients. Epilepsy lesions were on average ranked among the top 4 most suspicious detected clusters, thus demonstrating a promising performance for this very challenging detection task.

Figure 2 :

 2 Figure 2: General pipeline of our brain anomaly detection model. The training is shown in a yellow path, testing in a purple path. This pipeline illustrates the processing of monomodal T1-w MRI patches for clarity purpose. In the multimodal setting, patches of the T1-w and FLAIR imaging modalities extracted at the same spatial location in the brain are combined as channels at the input of the deep unsupervised representation model.

Figure 3 :

 3 Figure 3: Siamese neural network composed of stacked convolutional autoencoders as subnetworks. The input consists of a pair of patches of 2 different subjects centered at the same spatial localization in the brain. The middle-layer representation is denoted by z.

Figure 4 :

 4 Figure 4: An example of post-processing on a patient. The first column shows the original slice centered at the lesion, the second column corresponds to the normalized distance map Fp and the last three columns are obtained by thresholding Fp at three different values and identifying the connected components as detected clusters.

  Performance of our brain anomaly detection modelThe brain anomaly detection model of figure was trained on the training dataset of the 69 healthy subjects as described in section 4.1. The validation dataset was used to define convergence of the rSN model and fit the parameter set Θ of the loss term in eq. 1. The hyperparameter α that controls the tradeoff between the two terms of the rSN loss term was varied among the values : 0.25 and 0.5. Detection performance of the model was evaluated on the series of 21 patients based on the metric defined in section 4.4.

Figure 5 :

 5 Figure 5: Encoder E and generator G used in the convolutional autoencoder, Wasserstein autoencoder and regularized siamese network multichannel architectures.

Figure 6 :

 6 Figure 6: Left : Siamese neural network composed of stacked convolutional autoencoders as subnetworks (CAE). Center: Convolutional Autoencoder (CAE) composed of an encoder E and decoder D. Right: Wasserstein autoencoder (WAE) composed of an encoder E, a decoder G and an adversary discriminator D. For all three models, the encoder and decoder have the same architecture described on figure 5.

Figure 7 :

 7 Figure 7: Output detection maps for patients D + , N -, B -and E -respectively ( + stands for MRI-positive patients, -for MRI-negative patients). Top row: Transverse slices centered at the lesion locations (highlighted in purple). Bottom row: Maximum intensity projections (MIP) of the cluster maps overlaid on the MRI transverse slices. The maps show the top 1, top 1, top 1 and top 8 clusters, respectively, corresponding to the rank of the true positive detections. The color scale indicates the ranking of the detected clusters, ranging from light yellow (for the most suspicious ones) to dark brown (for the lowest ranked ones).

Figure 8 :

 8 Figure 8: Comparative fROC curves estimated on the 21 patients of the test dataset. x-axis: number of false positive detections per patient, y-axis: sensitivity. a) Performance of the proposed detection model based on T1/FLAIR rSN representation and comparison with similar architecture trained on a single imaging modality (T1 or FLAIR). T1, FLAIR and T1/FLAIR rSNs are shown for two values of α. b) Comparison of our model with alternate architecture based on an intermediate fusion strategy of T1 and FLAIR modalities using slimSimpleMKL. Five values of the parameter λ were considered for slimSimpleMKL. fROC curves are superimposed for λ ≥ 0.1. c) Comparison of our model with alternate architectures based on two different unsupervised representation models, CAE and WAE.

Table 2 :

 2 State-of-the art performances for the detection of FCD lesions in brain T1-w MRI. First column: Sensitivity and corresponding average FP rate inside parentheses. Note that the false positive rate is estimated on a normal control group. Second column: Size of the patient test dataset.

		Sensitivity (FP)	Nb of Patients
	Gill et al. (2018)	0.87-0.90 (4)	102 M RI +
	Gill et al. (2017)	0.83 (4)	41 M RI +
	Tan et al. (2018)	0.82 (26)	28 M RI +
	Jin et al. (2018)	0.74 (90%)	44 M RI + and 17 M RI -
	Jin et al. (2018)	0.53 (90%)	17 M RI -
	Hong et al. (2014b) 0.74 (95-100%)	17 M RI +

Table 3 :

 3 Performance of our brain anomaly detection model compared with state-of-the art performances in similar experimental condition for the detection of M RI -lesions. First column: Sensitivity and corresponding average false positive rate inside parentheses. Second column: Size of the patient test dataset. Third column: MRI sequences Gill et al. (2017); Jin et al. (2018); Tan et al. (2018)), thus restricting their analysis to the gray matter area. Recent studies by Ahmed et al. (2016) and El Azami et al. (2016) both restricted their evaluation to M RI -patients and use the same performance metrics as ours thus allowing a more straightforward comparison as reported in table 3. The system proposed in[START_REF] Ahmed | Decrypting "cryptogenic" epilepsy: Semi-supervised hierarchical conditional random fields for detecting cortical lesions in mri-negative patients[END_REF] based on SBM features coupled with semi-supervised hierarchical conditional random fields achieves 70% sensitivity among the top 10 detections per scan, while in El Azami et al. (2016), our model based on morphometric and intensity features coupled with a oc-SVM classifier allows achieving the same 70% sensitivity with an average of 4 false positives per scan.The slightly superior performance achieved by these two methods compared to our model is likely to be explained by the use of handcrafted features targeting FCD lesions. To emphasize this point, we implemented the method proposed in El Azami et al. (

		Sensitivity (FP) Nb of M RI -Patients	MRI
	Ahmed et al. (2016)	0.70 (9)	20	T1
	El Azami et al. (2016)	0.70 (4)	10	T1
	our implementation of El Azami et al. (2016)	0.52 (9)	18	T1
	our model	0.62 (9)	18	T1/FLAIR
	our model with T1 rSN	0.39 (9)	18	T1
	our model with FLAIR rSN	0.50 (9)	18	FLAIR
	Ahmed et al. (2016);

Hereafter, a rSN with α = m is a shorthand for αmax = m.

Acknowledgements

This work received funds from the French Foundation for Research in Epilepsy (FFRE) and the Region Auvergne-Rhône-Alpes through the TADALOT project. It was performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).

The authors sincerely thank Jiazheng Chai for his valuable contribution to the implementation of the WAE model. The authors would also like to thank the Institute of Informatics, Slovak Academy of Sciences for access to GPU resources through the VIP portal and technical support.

Comparison with unsupervised models based on handcrafted features

We compare the performance of our brain anomaly detection model to those obtained with two configurations based on handcrafted features.

The first configuration is a general linear model (GLM) (SPM analysis) learned on feature maps derived from T1w images using three settings -1. junction contrast, 2. extension contrast and 3. the conjunction of both contrastsfor a statistical threshold of 0.001 as done in El Azami et al. (2016). These feature maps model the junction between the gray and white matters as described in Huppertz et al. (2005) and Wagner et al. (2011). For a fair comparison, the same clustering and ranking procedures as performed for our anomaly detection model (see section 3.4) were applied and the top 10 most suspicious clusters were considered. The results are summarized in table 1. While extension contrast detects one additional lesion compared with our monomodal architecture based on T1w MRI (T1 rSN), the combination of junction and extension contrasts results in an inferior performance. This shows that it is not trivial to perform a multivariate analysis within this approach, retaining the best performance obtained in the univariate case. We should also note that without applying the ranking method, the original SPM implementation produces much more false positive detections without any significant change in sensitivity.

The second configuration whose performance is also summarized in table 1 is the model proposed in El Azami et al. (2016) where a oc-SVM is built per voxel using the same junction and extension maps, as above. This comparison Preprint submitted to Journal of L A T E X Templates November 24, 2020 reveals the advantage of the representations learnt with the rSN network versus the handcrafted features. As it can be seen, the maximum sensitivity achieved with the handcrafted features is 8/21 while 9/21 with our system on T1-w MRI, for 9 FPs. The gap becomes even more remarkable when the FLAIR modality is considered, leading our system to achieve 62% sensitivity.