Benoit Da Mota
email: benoit.damota@inria.fr

Radu Tudoran

Alexandru Costan

Ga Ël Varoquaux

Goetz Brasche

Patricia Conrod

Herve Lemaitre

Tomas Paus

Marcella Rietschel

Vincent Frouin

Jean-Baptiste Poline

Gabriel Antoniu

Bertrand Thirion
email: bertrand.thirion@inria.fr

Imagen Consortium

Frontiers in Neuroinformatics

Keywords: machine learning, neuroimaging-genetic, cloud computing, fMRI, heritability

Brain imaging is a natural intermediate phenotype to understand the link between genetic information and behavior or brain pathologies risk factors. Massive efforts have been made in the last few years to acquire high-dimensional neuroimaging and genetic data on large cohorts of subjects. The statistical analysis of such data is carried out with increasingly sophisticated techniques and represents a great computational challenge. Fortunately, increasing computational power in distributed architectures can be harnessed, if new neuroinformatics infrastructures are designed and training to use these new tools is provided. Combining a MapReduce framework (TomusBLOB) with machine learning algorithms (Scikitlearn library), we design a scalable analysis tool that can deal with non-parametric statistics on high-dimensional data. End-users describe the statistical procedure to perform and can then test the model on their own computers before running the very same code in the cloud at a larger scale. We illustrate the potential of our approach on real data with an experiment showing how the functional signal in subcortical brain regions can be significantly fit with genome-wide genotypes. This experiment demonstrates the scalability and the reliability of our framework in the cloud with a two weeks deployment on hundreds of virtual machines.

INTRODUCTION

Using genetics information in conjunction with brain imaging data is expected to significantly improve our understanding of both normal and pathological variability of brain organization. It should lead to the development of biomarkers and in the future personalized medicine. Among other important steps, this endeavor requires the development of adapted statistical methods to detect significant associations between the highly heterogeneous variables provided by genotyping and brain imaging, and the development of software components with which large-scale computation can be done.

In current settings, neuroimaging-genetic datasets consist of a set of i) genotyping measurements at given genetic loci, such as Single Nucleotide Polymorphisms (SNPs) that represent a large amount of the genetic between-subject variability, and ii) quantitative measurements at given locations (voxels) in threedimensional images, that represent e.g. either the amount of functional activation in response to a certain task or an anatomical feature, such as the density of grey matter in the corresponding brain region. These two sets of features are expected to reflect differences in brain organization that are related to genetic differences across individuals.

Most of the research efforts so far have been focused on designing association models, while the computational procedures used to run these models on actual architectures have not been considered carefully. Voxel intensity and cluster size methods have been used for genome-wide association studies (GWAS) [START_REF] Stein | Voxelwise genome-wide association study (vGWAS)[END_REF], but the multiple comparisons problem most often does not permit to find significant results, despite efforts to estimate the effective number of tests [START_REF] Gao | Avoiding the high Bonferroni penalty in genome-wide association studies[END_REF] or by paying the cost of a permutation test [START_REF] Da Mota | A fast computational framework for genome-wide association studies with neuroimaging data[END_REF]. Working at the genes level instead of SNPs [START_REF] Hibar | Voxelwise gene-wide association study (vGeneWAS): multivariate gene-based association testing in 731 elderly subjects[END_REF][START_REF] Ge | Increasing power for voxelwise genome-wide association studies: the random field theory, least square kernel machines and fast permutation procedures[END_REF]) is a promising approach, especially if we are looking at monogenic (or few causal genes) diseases.

For polygenic diseases, gains in sensitivity might be provided by multivariate models in which the joint variability of several genetic variables is considered simultaneously. Such models are thought to be more powerful [START_REF] Vounou | Discovering genetic associations with high-dimensional neuroimaging phenotypes: A sparse reduced-rank regression approach[END_REF][START_REF] Bunea | Penalized least squares regression methods and applications to neuroimaging[END_REF][START_REF] Kohannim | Boosting power to detect genetic associations in imaging using multi-locus, genome-wide scans and ridge regression[END_REF][START_REF] Meinshausen | Stability selection[END_REF][START_REF] Floch | Significant correlation between a set of genetic polymorphisms and a functional brain network revealed by feature selection and sparse partial least squares[END_REF], because they can express more complex relationships than simple pairwise association models. The cost of unitary fit is high due to high-dimensional, potentially nonsmooth optimization problems and various cross-validation loops needed to optimize the parameters; moreover, permutation testing is necessary to assess the statistical significance of the results of such procedures in the absence of analytical tests. Multivariate statistical methods require thus many efforts to be tractable for this problem on both the algorithmic and implementation side, including the design of adapted dimension reduction schemes. Working in a distributed context is necessary to deal efficiently with the memory and computational loads.

Today, researchers have access to many computing capabilities to perform data-intensive analysis. The cloud is increasingly used to run such scientific applications, as it offers a reliable, flexible, and easy to use processing pool [START_REF] Vaquero | A break in the clouds: towards a cloud definition[END_REF][START_REF] Juve | An evaluation of the cost and performance of scientific workflows on amazon ec2[END_REF][START_REF] Jackson | Seeking supernovae in the clouds: a performance study[END_REF][START_REF] Hiden | Developing cloud applications using the e-science central platform[END_REF]. The MapReduce paradigm [START_REF] Chu | Map-reduce for machine learning on multicore[END_REF][START_REF] Dean | MapReduce: simplified data processing on large clusters[END_REF] is the natural candidate for these applications, as it can easily scale the computation by applying in parallel an operation on the input data (map) and then combine these partials results (reduce). However, some substantial challenges still have to be addressed to fully exploit the power of cloud infrastructures, such as data access, as it is currently achieved through high latency protocols, which are used to access the cloud storage services (e.g. Windows Azure Blob). To sustain geographically distributed computation, the storage system needs to manage concurrency, data placement and inter-site data transfers.

We propose an efficient framework that can manage inferences on neuroimaging-genetic studies with several phenotypes and permutations. It combines a MapReduce framework (TomusBLOB, [START_REF] Costan | TomusBlobs: Scalable Data-intensive Processing on Azure Clouds[END_REF]) with machine learning algorithms (Scikit-learn library) to deliver a scalable analysis tool. The key idea is to provide end-users the capability to easily describe the statistical inference that they want to perform and then to test the model on their own computers before running the very same code in the cloud at a larger scale. We illustrate the potential of our approach on real data with an experiment showing how the functional signal in subcortical brain regions of interest (ROIs) can be significantly predicted with genome-wide genotypes. In section 2, we introduce methodological prerequisites, then we describe our generic distributed machine learning approach for neuroimaging-genetic investigations and we present the cloud infrastructure. In section 3, we provide the description of the experiment and the results of the statistical analysis.

MATERIALS AND METHODS

NEUROIMAGING-GENETIC STUDY

Neuroimaging-genetic studies test the effect of genetic variables on imaging target variables in presence of exogenous variables. The imaging target variables are activation images obtained through functional Magnetic Resonance Imaging (fMRI), that yield a standardized effect related to experimental stimulation at each brain location of a reference brain space. For a study involving n subjects, we generally consider the following model:

Y = Xβ 1 + Zβ 2 + ǫ,
where Y is a n×p matrix representing the signal of n subjects described each by p descriptors (e.g. voxels or ROIs of an fMRI contrast image), X is the n × q 1 set of q 1 explanatory variables and Z the n × q 2 set of q 2 covariates that explain some portion of the signal but are not to be tested for an effect. β 1 and β 2 are the fixed coefficients of the model to be estimated, and ǫ is some Gaussian noise. X contains genetic measurements and variables in Z can be of any type (genetic, artificial, behavioral, experimental, . . . ).

The standard approach. It consists in fitting p Ordinary Least Square (OLS) regressions, one for each column of Y, as a target variable, and each time perform a statistical test (e.g. F-test) and interpret the results in term of significance (p-value). This approach suffers from some limitations. First, due to a low signal-to-noise ratio and a huge number of tests, this approach is not sensitive. Moreover, the statistical score only reflects the univariate correlation between a target and a set of q 1 explanatory variables, it does not inform on their predictive power when considered jointly. Secondly, with neuroimaging data as a signal, we are not in a case vs. control study. It raises the question whether the variability in a population can be imputed to few rare genetic variants or if it is the addition of many small effects of common variants. Unfortunately, the model holds only if n ≫ (q 1 + q 2 ), which is not the case with genome-wide genotypes.

Heritability assessment. The goal of our analysis is to estimate the proportion of differences in a trait between individuals due to genetic variability. Heritability evaluation traditionally consists in studying and comparing homozygous and dizygous twins, but recently it has been shown that it can be estimated using genome-wide genotypes (Yang et al., 2011a;[START_REF] Lee | Estimating missing heritability for disease from genome-wide association studies[END_REF][START_REF] Lippert | Fast linear mixed models for genome-wide association studies[END_REF]. For instance, common variants are responsible of a large portion of the heritability of human height [START_REF] Yang | Common snps explain a large proportion of the heritability for human height[END_REF] or schizophrenia [START_REF] Lee | Estimating the proportion of variation in susceptibility to schizophrenia captured by common snps[END_REF]. These results show that the variance explained by each chromosome is proportional to its length. As we consider fMRI measurements in an unsupervised setting (no disease), this suggests to use regression models that do not enforce sparsity. Like the standard approach, heritability has some limitations. In particular, the estimation of heritability requires large sample sizes to have an acceptable standard error (at least 4000 according to [START_REF] Lee | Estimating the proportion of variation in susceptibility to schizophrenia captured by common snps[END_REF]). Secondly, the heritability is the ratio between the variance of the trait and the genetic variance in a population. Therefore, for a given individual, a trait with an heritability at 0.6 does not mean it can be predicted at 60% on average with the genotype. It means that a fraction of the phenotype variability is simply explained by the average genetic structure of the population of interest.

High-dimensional statistics. The key point of our approach is to fit a model on training data (train set) and evaluate its goodness on unseen data (test set). To stabilize the impact of the sets for training and testing, a cross-validation loop is performed, yielding an average prediction score over the folds. This score yields a statistic value and a permutation test is performed to tabulate the distribution of this statistic under the null hypothesis and to estimate its significance (p-value). In practice, this corresponds to swapping the labels of the observations. As a prediction metric we generally choose the coefficient of determination (R 2 ), which is the ratio between the variance of the prediction and the variance of the phenotypes in the test set. If we consider all the genotypes at the same time, this approach is clearly related to heritability, but focuses on the predictive power of the model and its significance. Through cross-validation, the estimation of the CV -R 2 with an acceptable standard error does not require as large sample sizes as for the estimation of heritability (Yang et al., 2011b).

CV -R 2 = 1 -mean (train,test)∈split(n) Y test -X test β train 1 -Z test β train 2 2 Y test -Z test β train 2 2

GENERIC PROCEDURE FOR DISTRIBUTED MACHINE LEARNING

If one just wants to compute the prediction score for few phenotypes, a multicore machine should be enough. But, if one is interested in the significance of this prediction score, one will probably need a computers farm (cloud, HPC cluster, etc.) Our approach consists in unifying the description and the computation for neuroimaging-genetic studies to scale from the desktop computer to the supercomputing facilities. The description of the statistical inference is provided by a descriptive configuration in humanreadable and standard format: JSON (JavaScript Object Notation). This format requires no programing skills and is far easier to process as compared to the XML (eXtensible Markup Language) format. In a sense, our approach extends the Scikit-learn library (cf. next paragraph) for distributed computing, but focuses on a certain kind of inferences for neuroimaging-genetic studies. The next paragraphs describe the strategy, framework and implementation used to meet the heritability assessment objective. Scikit-learn is a popular machine learning library in Python [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] designed for a multicore station. In the Scikit-learn vocabulary, an estimator is an object that implements a fit and a predict method. For instance a Ridge object (lines 12-13 of Figure 1) is an estimator that computes the coefficients of the ridge regression model on the train set and uses these coefficients to predict data from the test set. If this object has a transform method, it is called a transformer. For instance a SelectKbest object (lines 10-11 of Figure 1) is a transformer that modifies the input data (the design matrix X) by returning the K best explanatory variables w.r.t. a scoring function. Scikit-learn defines a Pipeline (lines 8-13 of Figure 1) as the combination of several transformers and an final estimator: It creates a combined estimator. Model selection procedures are provided to evaluate with a cross-validation the performance of an estimator (eg. cross val score) or to select parameters on a grid (eg. GridSearchCV).

Permutations and covariates. Standard machine learning procedures have not been designed to deal with covariates (such as those assembled in the matrix Z), which have to be considered carefully in a permutation test [START_REF] Anderson | Permutation tests for linear models[END_REF]. For the original data, we fit an Ordinary Least Square (OLS) model between Y and Z, then we consider the residuals of the regression (denoted R Y |Z ) as the target for the machine learning estimator. For the permutation test, we permute R Y |Z (the permuted version is denoted R Y |Z * ), then we fit an OLS model between R Y |Z * and Z, and we consider the residuals as the target for the estimator [START_REF] Anderson | Permutation tests for linear models[END_REF]. The goal of the second OLS on the permuted residuals is to provide an optimal approximation (in terms of bias and computation) of the exact permutation tests while working on the reduced model. Generic problem. We identify a scheme common to the different kinds of inference that we would like to perform. For each target phenotype we want to compute a prediction score in the presence of covariates or not and to evaluate its significance with a permutation test. Scikit-learn algorithms are able to execute on multiple CPU cores, notably cross-validation loop, so a task will be executed on a multicore machine: cluster nodes or multicore virtual machine (VM). As the computational burden of different machine learning algorithms is highly variable, owing to the number of samples and the dimensionality of the data, we thus have to tune the number of tasks and their average computation time. An optimal way to tune the amount of work is to perform several permutations on the same data in a given task to avoid I/O bottlenecks. Finally, we put some constraints on the description of the machine learning estimator and the cross validation scheme:

• The prediction score is computed using the Scikit-learn cross val score function and the folds for this cross validation loop are generated with a ShuffleSplit object. • An estimator is described with a Scikit-learn pipeline with one or more steps. prediction score> 1 {"extract_cov": true, 2 "n_perm_total": 10000, 3 "n_perm_per_mapper": 1, 4 "cross_val_score": { 5 "score_func": "DYNAMIC_IMPORT::sklearn.metrics.r2_score"}, 6 "ShuffleSplit": { 7 "test_size": 0.2, "random_state": 0, "n_iter": 10}, 8 "pipeline": [ 9 ["FastFilterColinear", "gstat.data.utils.FastFilterColinear", {}], 10 ["SelectKBest", "sklearn.feature_selection.SelectKBest", { 11 "score_func": "DYNAMIC_IMPORT::sklearn.feature_selection.f_regression"}],

12 ["Ridge", "sklearn.linear_model.Ridge", { 13 "fit_intercept": true}]],

14 "GridSearchCV": ["sklearn.grid_search.GridSearchCV", {}, [{ 15 "SelectKBest__k" : [10, 100, 1000],

16 "Ridge__alpha" : [0.0001, 0.001, 0.01, 0.1, 1.]}]]} Figure 1. (Top) representation of the computational framework: given the data, a permutation and a phenotype index together with a configuration file, a set of computations are performed, that involve two layers of cross-validation for setting the hyper-parameters and evaluate the accuracy of the model. This yields a statistical score associated with the given phenotype and permutation. (Bottom) Example of complex configuration file that describes this set of operations. General parameters (Lines 1-3): The model contains covariates, the permutation test makes 10,000 iterations and only one permutation is performed in a task. Prediction score (Lines 4-7): The metrics for the cross-validated prediction score is R 2 , the cross-validation loop makes 10 iterations, 20% of the data are left out for the test set and the seed of the random generator was set to 0. Estimator pipeline (Lines 8-13): The first step consists in filtering collinear vectors, the second step selects the K best features and the final step is a ridge estimator. Parameters selection (Lines 14-16): 2 parameters of the estimator have to be set: the K for the SelectKBest and the alpha of the Ridge regression. A set of 3 × 5 parameters are evaluated.

• Python can dynamically load modules such that a program can execute functions that are passed in a string or a configuration file. To notify that a string contains a Python module and an object or function to load, we introduce the prefix DYNAMIC IMPORT:: • To select the best set of parameters for an estimator, model selection is performed using Scikit-learn

GridSearchCV and a 5-folds inner cross-validation loop.

Full example (cf. script in Figure 1):

• General parameters (Lines 1-3): The model contains covariates, the permutation test makes 10,000 iterations and only one permutation is performed in a task. 10,000 tasks per brain target phenotypes will be generated. • Prediction score (Lines 4-7): The metrics for the cross-validated prediction score is R 2 , the crossvalidation loop makes 10 iterations, 20% of the data are left out for the test set and the seed of the random generator was set to 0. • Estimator pipeline (Lines 8-13): The first step consist in filtering collinear vectors, the second step selects the K best features and the final step is a ridge estimator. • Parameters selection (Lines 14-16): 2 parameters of the estimator have to be set: the K for the SelectKBest and the alpha of the Ridge regression. A set of 3 × 5 parameters are evaluated.

User Browser

Management Deployment

Global Reducer Splitter

THE CLOUD COMPUTING ENVIRONMENT

Although researchers have relied mostly on their own clusters or grids, clouds are raising an increasing interest [START_REF] Juve | An evaluation of the cost and performance of scientific workflows on amazon ec2[END_REF][START_REF] Jackson | Seeking supernovae in the clouds: a performance study[END_REF][START_REF] Ghoshal | I/o performance of virtualized cloud environments[END_REF][START_REF] Simmhan | Bridging the gap between desktop and the cloud for escience applications[END_REF][START_REF] Hiden | Developing cloud applications using the e-science central platform[END_REF]. While shared clusters or grids often imply a quota-based usage of the resources, those from clouds are owned until they are explicitly released by the user. Clouds are easier to use since most of the details are hidden to the end user (eg. network physical implementation). Depending on the characteristics of the targeted problem, this is not always an advantage (eg. collective communications). Last but not least, clouds avoid owning expensive infrastructures -and associated high cost for buying and operating-that require technical expertise. The cloud infrastructure is composed of multiple data centers, which integrate heterogeneous resources that are exploited seamlessly. For instance, the Windows Azure cloud has 5 sites in United States, 2 in Europe and 3 in Asia. As resources are granted on-demand, the cloud gives the illusion of infinite resources. Nevertheless, cloud data centers face the same load problems (e.g. workload balancing, resource idleness, etc.) as traditional grids or clusters.

In addition to the computation capacity, clouds often provide data-related services, like object storage for large datasets (e.g. S3 from Amazon or Windows Azure Blob) and queues for short message communication.

NEUROIMAGING-GENETICS COMPUTATION IN THE CLOUD

In practice, the workload of the A-Brain application1 is more resource demanding than the typical cloud applications and could induce two undesirable situations: 1) other clients do not have enough resource to lease on-demand in a particular data center; 2) the computation creates performance degradations for other applications in the data center (e.g. by occupying the network bandwidth, or by creating high number of concurrent requests on the cloud storage service). Therefore, we divide the workload into smaller sub-problems and we select the different datacenters in collaboration with the cloud provider.

For balancing the load of the A-Brain application, the computation was distributed across 4 deployments in the two biggest Windows Azure datacenters. In the cloud context, a deployment denotes a set of leased resources, which are presented to the user as a set of uniform machines, called compute nodes. Each deployment is independent and isolated from the other deployments. When a compute node starts, the user application is automatically uploaded and executed. The compute nodes of a deployment belong to the same virtual private network and communicate with the outside world or other deployments either through public endpoints or using the cloud storage services (i.e. Windows Azure Blob or Queue).

TomusBlobs [START_REF] Costan | TomusBlobs: Scalable Data-intensive Processing on Azure Clouds[END_REF] is a data management system designed for concurrency-optimized PaaS-level (Platform as a Service) cloud data management. The system relies on the available local storage of the compute nodes in order to share input files and save output files. We built a processing framework (called TomusMapReduce) derived from MapReduce [START_REF] Chu | Map-reduce for machine learning on multicore[END_REF][START_REF] Dean | MapReduce: simplified data processing on large clusters[END_REF] on top of TomusBlobs, such that it leverages its benefits by collocating data with computation. Additionally, the framework is restricted to associative and commutative reduction procedures (Map-IterativeReduce model [START_REF] Tudoran | Mapiterativereduce: a framework for reduction-intensive data processing on azure clouds[END_REF]) in order to allow efficient out-of-order and parallel processing for the reduce phase. Although MapReduce is designed for single cluster processing, the latter constraint enables straightforward geographically distributed processing. The hierarchical MapReduce (which is described in [START_REF] Costan | TomusBlobs: Scalable Data-intensive Processing on Azure Clouds[END_REF]) aggregates several deployments with MapReduce engines and a last deployment that contains a MetaReducer, that computes the final result, and a Splitter, that partitions the data and manages the overall workload in order to leverage data locality. Job descriptions are sent to the MapReduce engines via Windows Azure Queue and the MetaReducer collects intermediate results via Windows Azure Blob. For our application, we use the Windows Azure Blob storage service instead of TomusBlobs for several reasons: 1) concurrency-optimized capabilities are not relevant here; 2) for a very long run, it is better to rely on a proven storage; 3) TomusBlob storage does not support yet multi-deployments setting. An overview of the framework is shown in Figure 2.

For our application, the Map step yields a prediction score for an image phenotype and a permutation, while the reduce step consists in collecting all results to compute statistic distribution and corrected p-values. The reduce operation is trivially commutative and associative as it consists in searching the maximum of the statistic for each permutation [START_REF] Westfall | Resampling-based multiple testing : examples and methods for P-value adjustment[END_REF]. The upper part of Figure 1 gives an overview of the generic mapper.

IMAGEN: A NEUROIMAGING-GENETIC DATASET

IMAGEN is a European multi-centric study involving adolescents [START_REF] Schumann | The imagen study: reinforcement-related behaviour in normal brain function and psychopathology[END_REF]. It contains a large functional neuroimaging database with fMRI associated with 99 different contrast images for 4 protocols in more than 2,000 subjects, who gave informed signed consent. Regarding the functional neuroimaging data, we use the Stop Signal Task protocol [START_REF] Logan | On the ability to inhibit thought and action: A users' guide to the stop signal paradigm[END_REF] (SST), with the activation during a [go wrong] event, i.e. when the subject pushes the wrong button. Such an experimental contrast is likely to show complex mental processes (inhibition failure, post-hoc emotional reaction of the subject), that may be hard to disentangle. Our expectation is that the amount of Blood Oxygen-Level Dependent (BOLD) response associated with such events provides a set of global markers that may reveal some heritable psychological traits of the participants. Eight different 3T scanners from multiple manufacturers (GE, Siemens, Philips) were used to acquire the data. Standard preprocessing, including slice timing correction, spike and motion correction, temporal detrending (functional data) and spatial normalization (anatomical and functional data), were performed using the SPM8 software and its default parameters; functional images were resampled at 3mm resolution. All images were warped in the MNI152 coordinate space. Obvious outliers detected using simple rules such as large registration or segmentation errors or very large motion parameters were removed after this step. BOLD time series was recorded using Echo-Planar Imaging, with TR = 2200 ms, TE = 30 ms, flip angle = 75 • and spatial resolution 3mm × 3mm × 3mm. Gaussian smoothing at 5mm-FWHM was finally added. Contrasts were obtained using a standard linear model, based on the convolution of the time course of the experimental conditions with the canonical hemodynamic response function, together with standard high-pass filtering procedure and temporally auto-regressive noise model. The estimation of the first-level was carried out using the SPM8 software. T1-weighted MPRAGE anatomical images were acquired with spatial resolution 1mm × 1mm × 1mm, and gray matter probability maps were available for 1986 subjects as outputs of the SPM8 New Segmentation algorithm applied to the anatomical images. A mask of the gray matter was built by averaging and thresholding the individual gray matter probability maps. More details about data preprocessing can be found in [START_REF] Thyreau | Very large fMRI study using the IMAGEN database: sensitivity-specificity and population effect modeling in relation to the underlying anatomy[END_REF]. 1 {"extract_cov": true, 2 "n_perm_total": 10000, 3 "n_perm_per_mapper": 5, 4 "cross_val_score": { 5 "score_func": "DYNAMIC_IMPORT::sklearn.metrics.r2_score"}, 6 "ShuffleSplit": { 7 "test_size": 0.2, "random_state": 0, "n_iter": 10}, 8 "pipeline": [ 9 ["SelectKBest", "sklearn.feature_selection.SelectKBest", { 10 "score_func": "DYNAMIC_IMPORT::gstat.stats.utils.f_regression", 11 "k": 50000}],

12 ["Ridge", "sklearn.linear_model.Ridge", { 13 "fit_intercept": true, "alpha": 0.0001}]]} Figure 3. Configuration used for the experiment. (Lines 1-3): covariates, 10,000 permutations and 5 permutations per computation unit (mapper). (Lines 4-7): 10-folds cross-validated R 2 .(Lines 9-11): The first step of the pipeline is an univariate features selection (K=50,000). This step is used as a dimension reduction so that the next step fits in memory. (Lines 12-13): The second and last step is the ridge estimator with a low penalty (alpha=0.0001).

DNA was extracted from blood samples using semi-automated process. Genotyping was performed genome-wide using Illumina Quad 610 and 660 chips, yielding approximately 600,000 autosomic SNPs. 477,215 SNPs are common to the two chips and pass plink standard parameters (Minor Allele Frequency > 0.05, Hardy-Weinberg Equilibrium P < 0.001, missing rate per SNP < 0.05).

AN APPLICATION AND RESULTS

THE EXPERIMENT

The aim of this experiment is to show that our framework has the potential to explore links between neuroimaging and genetics. We consider an fMRI contrast corresponding to events where subjects make motor response errors ([go wrong] fMRI contrast from a Stop Task Signal protocol). Subjects with too many missing voxels or with bad task performance were discarded. Regarding genetic variants, 477,215 SNPs were available. Age, sex, handedness and acquisition center were included in the model as confounding variables. Remaining missing data were replaced by the median over the subjects for the corresponding variables. After applying all exclusion criteria 1,459 subjects remained for analysis. Analyzing the whole brain with all the genetic variants remains intractable due to the time and memory requirements and dimension reduction techniques have to be employed. Prior neuroimaging dimension reduction. In functional neuroimaging, brain atlases are mainly used to provide a low-dimensional representation of the data by considering signal averages within groups of neighboring voxels. In this experiment we focus on the subcortical nuclei using the Harvard-Oxford subcortical atlas. We extract the functional signal of 14 regions of interest, 7 in each hemisphere: thalamus, caudate, putamen, pallidum, hippocampus, amygdala and accumbens (see Figure 4). White matter, brain stem and ventricles are of no interest for functional activation signal and were discarded. This prior dimension reduction decreases the number of phenotypes from more than 50,000 voxels to 14 ROIs.

Configuration used (cf. script in Figure 3):

• (Lines 1-3): covariates, 10,000 permutations and 5 permutations per computation unit (mapper).

• (Lines 4-7): 10-folds cross-validated R 2 .

• (Lines 9-11): The first step of the pipeline is an univariate features selection (K=50,000). This step is used as a dimension reduction so that the next step fits in memory. • (Lines 12-13): The second and last step is the ridge estimator with a low penalty (alpha=0.0001).

The goal of the experiment described by this configuration file is to evaluate how the 50,000 mostly correlated genetic variants, once taken together, are predictive of each ROI and to associate a p-value with these prediction scores. Note that more than 50,000 covariates would not fit into memory. This configuration generates 28,000 map tasks (14 × 10000/5), but we can set to 1 the number of permutations per task, which means that the computation can use up to 140,000 multicore computers in parallel, and thus millions of CPU cores. The cloud experimental setup. The experiment was performed using the Microsoft Windows Azure PaaS cloud in the North and West US datacenters, that were recommended by the Microsoft team for their capacity. We use the Windows Azure storage services (Blob and Queue) in both datacenters in order to take advantage of the data locality. Due to our memory requirements, the Large VM type (4 CPU cores, 7GB of memory and 1,000GB of disk) is the best fit regarding the Azure VMs offer2 .

TomusBlobs. We set up 2 deployments in each of the 2 recommended sites for a total of 4 deployments. It used 250 large VM nodes, totalizing 1,000 CPUs: each of the 3 MapReduce engines deployments had 82 nodes and the last deployment used 4 nodes. The reduction process was distributed in approximately 600 reduce jobs.

RESULTS

Cloud aspects. The experiment timespan was 14 days. The processing time for a single map job is approximately 2 hours. There are no noticeable time differences between the execution times of the map jobs with respect to the geographical location. In large infrastructures like the clouds, failures are possible and applications need to cope with this. In fact, during the experiment the Azure services became temporary inaccessible3 , due to a failure of a secured certificate. Despite this problem, the framework was able to handle the failure with a fault tolerance mechanism which suspended the computation until all Azure services became available again. The monitoring mechanism of the Splitter, that supervises the computation progress, was able to restore aborted jobs. The IterativeReduce approach eliminates the implicit barrier between mappers and reducers, but yields negligible gains due to the huge workload of the mappers. The effective cost of the experiment was approximately equal to 210,000 hours of sequential computation, which corresponds to almost 20,000$ (VM pricing, storage and outbound traffic).

Application side. Figure 4 shows a summary of the results. Despite the fact that some prediction scores are negative, the activation signal in each ROI is fit significantly better than chance using the 50,000 best genetic variants over the 477,215. The mean BOLD signal is better predicted in the left and right thalamus. The distribution of the CV -R 2 is also very informative, showing that by chance the mean prediction score is negative (familywise-error corrected or not). While this phenomenon is somewhat counter-intuitive within the framework of classical statistics, it should be pointed out that the cross-validation procedure used here opens the possibility of negative R 2 : this quantity is by definition a model comparison statistic that takes the difference between a regression model with a non-informative model; in high-dimensional settings, a poorly fitting linear model performs (much) worse than a non-informative model. Hence a model performing at chance gets a negative score: This is actually what happens systematically when the association between y and X is broken by the permutation procedure, even if we consider the supremum over many statistical tests [START_REF] Westfall | Resampling-based multiple testing : examples and methods for P-value adjustment[END_REF]. A slightly negative value can thus be the marker of a significant association between the variables of interest. Twin and SNP-based studies suggest high heritability of structural brain measures, such as total amount of gray and white matter, overall brain volume and addiction-relevant subcortical regions. Heritability estimates for brain measures are as high as 0.89 [START_REF] Kremen | Genetic and environmental influences on the size of specific brain regions in midlife: the vetsa mri study[END_REF] or even up to 0.96 [START_REF] Van Soelen | Brain scale: brain structure and cognition: an adolescent longitudinal twin study into the genetic etiology of individual differences[END_REF] and subcortical regions appear to be moderately to highly heritable. One recent study on subcortical volumes [START_REF] Den Braber | Heritability of subcortical brain measures: A perspective for future genome-wide association studies[END_REF] reports highest heritability estimates for the thalamus (0.80) and caudate nucleus (0.88) and lowest for the left nucleus accumbens (0.44). Despite the fact that the CV -R 2 metric is not exactly an heritability measurement, our metric evaluates the predictability of the fitted model (i.e. how well it predicts the activation signal of a brain region with genetic measurements on unseen data) which is a good proxy for heritability. Thus, our results confirm that brain activation signals are an heritable feature in subcortical regions. These experiments can be used as a basis to further localize the genetic regions (pathways or genes) that are actually predictive of the functional activation. An important extension of the present work is clearly to extend this analysis to the cortical regions.

CONCLUSION

The quantitative evaluation of statistical models with machine learning techniques represents an important step in the comprehension of the associations between brain image phenotypes and genetic data. Such approaches require cross validation loops to set the hyper-parameters and to evaluate performances. Permutations have to be used to assess the statistical significance of the results, thus yielding prohibitively expensive analyses. In this paper, we present a framework that can deal with such a computational burden. It relies on two key points: i) it wraps the Scikit-learn library to enable coarse grain distributed computation. Yet it enforces some restrictions, i.e. it solves only a given class of problems (pipeline structure, cross-validation procedure and permutation test). The result is a simple generic code (few lines) that provides the user a quick way to conduct early, small-scale investigations on its own computer or at a larger scale on a high-performance computing cluster. With JSON we provide a standard format for the description of statistical inference so that no programming skills are required and so that it can be easily generated from a webpage form. ii) TomusBLOB permits to execute seamlessly the very same code on the Windows Azure cloud. We could also disable some parts of TomusBLOB to achieve a good compromise between the capabilities and the robustness. We demonstrate the scalability and the efficiency of our framework with a two weeks geographically distributed execution on hundreds of virtual machines.

The results confirm that brain activation signals are an heritable feature.

  Score (outer CV loop) <(Phenotype ID, Permutation ID), data> + configuration file <(Phenotype ID, Permutation ID),

Figure 2 .

 2 Figure 2. Overview of the multi site deployment of a hierarchical Tomus-MapReduce compute engine. 1) The end-user uploads the data and configures the statistical inference procedure on a webpage. 2) The Splitter partitions the data and manages the workload. The compute engines retrieves job information trough the Windows Azure Queues. 3) Compute engines perform the map and reduce jobs. The managment deployment is informed of the progression via the Windows Azure Queues system and thus can manage the execution of the global reducer. 4) The user downloads the results of the computation on the webpage of the experiment.

Figure 4 .

 4 Figure 4. Results of the real data analysis procedure. (Left) predictive accuracy of the model measured by cross-validation, in the 14 regions of interest, and associated statistical significance obtained in the permutation test. (Up right) distribution of the CV -R 2 at chance level, obtained through a permutation procedure. The distribution of the max over all ROIs is used to obtain the family-wise error corrected significance of the test. (Bottom right) outline of the chosen ROIs.

Frontiers in Neuroinformatics

http://www.irisa.fr/kerdata/abrain/Frontiers in Neuroinformatics

http://msdn.microsoft.com/fr-fr/library/windowsazure/dn197896.aspx

Azure Failure Incident: http://readwr.it/tAqFrontiers in Neuroinformatics

ACKNOWLEDGEMENT

This work was supported primarily by the Microsoft INRIA joint centre grant A-brain and secondarily by the Digiteo ICoGeN grant and the ANR grant ANR-10-BLAN-0128. HPC investigations were carried out using the computing facilities of the CEA-DSV and CATI cluster. The data were acquired within the IMAGEN project, which receives research funding from the E.U. Community's FP6, LSHM-CT-2007-037286. This manuscript reflects only the author's views and the Community is not liable for any use that may be made of the information contained therein.


Régis Mélin

Denis Feinberg
email: denis.feinberg@neel.cnrs.fr

Hervé Courtois

Ciprian Padurariu

Andreas Pfeffer

Jean Eudes Duvauchelle

Francois Lefloch

Thibaut Jonckheere

Jérôme Rech

Thierry Martin

B Douçot

D.C. Josephson transport by quartets and other Andreev resonances in superconducting bijunctions

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Josephson junctions are weak links connecting two superconductors [START_REF] Tinkham | Introduction to Superconductivity[END_REF]. Cooper pairs crossing between the two condensates establish phase coherence [START_REF] Josephson | [END_REF]. The transport properties depend on the phase difference ϕ between the two reservoirs. Zero voltage transport involves the d.c. Josephson current made of a steady flow of Cooper pairs, while at nonzero voltage, due to the relation dϕ dt = 2eV h , pair transport manifests through the a.c. Josephson effect with oscillating Cooper pairs. In transparent enough SN S junctions, transport occurs through Andreev reflections at the SN interfaces, forming Andreev electron-hole bound states at zero voltage. At nonzero voltage below the superconducting gap ∆, quasiparticle transport can occur through multiple Andreev reflections (MAR) [3]. Thus, dissipationless d.c. transport at V = 0 and dissipative pair-assisted quasiparticle transport at V = 0 are well-separated phenomena.

Recently, set-ups involving three superconductors coupled by a common link have been proposed theoretically [4,5,6,7,8,9] and recently achieved experimentally [10,11,12]. Transport is then governed by two (instead of one) phase or voltage variables. Coherent mechanisms coupling the three superconductors altogether become possible and provide new d.c. channels which involve correlated motion of Cooper pairs. The lowest order processes correspond to quartets, e.g, two pairs exiting simultaneously from a reservoir and crossing the bijunction, one towards each of the other reservoirs [8,9]. This process is deeply related to so-called Cooper pair splitting that occurs in an hybrid N SN bijunction. Cooper pairs are split from S by socalled Crossed Andreev Reflection, and two energy-and spin-entangled electrons are sent in the two separated normal metal reservoirs [13,14,15,16,17,18,19,20,21,22,23,24,25]. When those reservoirs are instead superconducting, such a process should be duplicated in order to produce two quasiparticles in each of the side superconductors, that recombine into Cooper pairs [8,9,12]. The basic process thus involves four fermions, and it originates from splitting two Cooper pairs from a given superconductor. It involves four (instead of two) Andreev reflections within the weak link N , and it exchanges opposite-spin electrons of the two split Cooper pairs. Here we review the present theory of quartet and multipair coherent transport, and summarize a recent experiment where transport signatures of quartets have been obtained in a triterminal Al -Cu bijunction [12]. Section 2 recalls the adiabatic argument for multipair transport. Section 3 presents theoretical results for a double dot bijunction. Section 4 summarizes the recent experiment and the quartet interpretation.

The adiabatic limit of multipair transport

Let us consider three superconductors S 0,a,b with voltages V 0 = 0, V a and V b , connected by a short metallic link or a quantum (double) dot (Figure 1). Defining the phases ϕ 0 = 0, ϕ a , ϕ b , the total free energy at equilibrium (V a,b = 0) can be written as a doubly periodic function F (ϕ a , ϕ b ). Fourier transforming and expressing the current in one lead (say, S a ), one obtains:

I a = 2e h ∂F ∂ϕ a = n,m I a (n, m) sin(nϕ a + mϕ b ) (1)
At nonzero voltages, one writes:

ϕ a = ϕ 0a + 2e h V a t , ϕ b = ϕ 0b + 2e h V b t (2)
where the phases at the origin of time should not be omitted. Within the adiabatic approximation, one plugs the time-dependent phases into the I a (ϕ a , ϕ b ) relation [START_REF] Tinkham | Introduction to Superconductivity[END_REF]. Thus, if V a,b satisfy nV a + mV b = 0, only the harmonics sin p(nϕ a + mϕ b ) of I a have a nonzero static value, and a d.c. current flows in S a , as well as in S 0 and S b . This current involves n (m) pairs crossing altogether to S a (S b ), and n + m pairs coming from S 0 . The simplest of these multipair processes involves quartets, e.g. n = m = 1. It is important to notice that the lowest-order Josephson processes taking one single pair from one to the other superconductor average to zero. From Equations (1,2), one finds that a d.c. multipair current in S a appears:

I d.c. a (ϕ 0a , ϕ 0b ) = p I a,max (p) sin p(nϕ 0a + mϕ 0b ) (3)
This shows the relevance of the phase combination (nϕ 0a + mϕ 0b ) which controls coherent d.c. transport by multipair processes. This phase can be emulated either by imposing suitable fluxes in a circuit containing the bijunction, or by sending a d.c. current in the bijunction. The first method is demanding as it implies the control of two phase variables. The second one is somewhat similar to driving a current into a Josephson junction and exploring its zero-voltage branch. An equivalent in the bijunction of this zero-voltage branch is a nV a + mV b = 0 line in the (V a , V b ) plane.

The above simple calculation is valid only for tunneling barriers and voltages well below the gap. Yet, in this case, the quartet process which is a second order process might be weak. On the contrary, if the bijunction is metallic and has a high transparency, quartet transport can be strong and even resonant. In this case, a full nonequilibrium calculation should be performed. In what follows, we explore two situations, first, a bijunction made of a double quantum dot, achievable with a carbon nanotoube or a nanowire, and second, with disordered long metallic junctions.

Quartet and multipair transport in a double quantum dot bijunction

Recent progress on Cooper pair splitting has arisen with the fabrication of hybrid N a D a S 0 D b N b bijunctions where the dots D are created within a single wall carbon nanotube [23] or a semiconducting nanowire [24,25]. Evidence for Cooper pair splitting comes out from measurement of the nonlocal conductance. Consider now a bijunction S a D a S 0 D b S b with a similar geometry, but with all leads being superconducting. For sake of simplicity, one considers one single electronic level ε a,b in each dot. The coupling between the dots and the superconductors is supposed to be symmetric and described by the level broadening Γ = πν(0)t 2 where ν is the normal metal density of states in the superconductors and t is a dot-superconductor hopping parameter. The adiabatic argument of Section 2 suggests that the multipair current depends on the phase combination nϕ 0a + mϕ 0b . The calculation of the current-phase characteristics is performed using nonequilibrium Keldysh Green's functions [9]. This yields the solution as an effective partition function for the dot "dressed" by tunneling into the superconductors. This approach is nonperturbative and provides exact results. Details about the formalism can be found in Ref. [26]. The parameters of the calculation (in units of the gap) are : the dot energies ε a,b , the coupling Γ, and the voltage V . Notice that the calculation gives the full current, which for a given commensurate combination of voltages comprises not only the coherent d.c. multipair current, but also MAR quasiparticle currents.

Open dot regime

One first considers a regime in which each dot mimics a metallic junction, achieved by placing energy levels out of resonance |ε a,b | > ∆, with large couplings Γ > ∆. The results for the largest resonances are represented on Figure 2 and compared to the phase-independent part of the quasiparticle current. One sees that the multipair features appear as strong resonances. The quartet dispersion has a negative slope at the origin and is thus a π-junction in terms of the "quartet phase" ϕ a + ϕ b . This is due to the exchange process occurring as two split Cooper pairs are emitted by S 0 and recombine into one pair in S a and one in S b .

Closed and nonresonant dot regime

Consider now small coupling and nonresonant dots, e.g. Γ < |ε a,b | < ∆. Figure 3 shows the current-phase profiles for various values of V , increasing from 0.09∆ to 0.8∆. The maximum current increases by several orders of magnitude and is maximum when V ∼ |ε a,b |. Moreover, above the maximum, currents I a and I b start to depart from each other, signalling the onset of a strong MAR quasiparticle current. Remarkably, this MAR current is phase-dependent, a specific feature of the resonant quartet behaviour, which can be ascribed to interferences between various amplitudes within the bijunction (Figure 4). This calculation validates the existence of quartet and higher-order resonances featuring a d.c. current of correlated Cooper pairs. These currents are functions of a well-defined combination of the phases, that becomes a relevant and observable variable. Phase coherence is thus established between the three superconductors, despite the presence of voltage biases. One should underline that instead of Andreev bound states forming in a junction at equilibrium, here the d.c. Josephson-like multipair resonances are due to Andreev resonances, and they coexist with dissipative MAR channels. This coexistence is a novel feature in superconductivity. As a corollary, the maximum "critical" quartet current depends on the applied voltage, in a dramatic way in the example shown above (Section 3.2).

Quartet anomalies in metallic bijunctions

Copper and Aluminum bijunctions have been fabricated [12] by shadow mask evaporation technique [27], leading to highly transparent and uniform SN interfaces.

They form three-terminal junctions with a T-shape normal conductor connecting three superconducting electrodes S 0 , S a and S b (Figure 5). The width and length of the normal metal are about 0.6 µm and 1 µm. Using a diffusion constant for copper D = 100 cm 2 /s, one gets a Thouless energy E T h = hD/L 2 = 7 µeV . The superconducting Al energy gap is ∆ = 170 µeV . These Josephson junctions are thus in the long junction regime defined as the Thouless energy E T h being much smaller than the gap ∆, or equivalently a junction length L larger than the superconductor coherence length ξ s = hD/∆. The diffusion time is τ D = L 2 /D ≃ 0.1ns is much smaller than the inelastic time τ in ≃ 1 ns at 100 mK. The SEM image shows a bijunction sample with a T-shape geometry.

Three-terminal differential resistances were measured using a specific experimental set-up including 3 SQUIDs as current amplifiers (Figure 5) [28]. The differential resistance ∂Va ∂Ia is plotted as a 2D map in the (V a , V b ) plane (Figure 6). First, it displays anomalies along the axis, that correspond to the "direct" Josephson effect between S 0 and S a or S b (the one between S a and S b is not visible due to the location of the ac lock-in). Second, clear features appear at voltage values such as V i + V j -2V k = 0, where (i, j, k) = (a, b, 0), (b, 0, a), (0, a, b), corresponding to quartets emitted respectively by S 0 , S a , S b . The transversal I(V ) traces shown in Figure 7 display comparable shape and amplitude with respect to the Josephson features obtained at V a = 0 or V b = 0.

An explanation in terms of the synchronization of ac Josephson oscillations in junctions S 0 S a and S 0 S b [29,30] is unlikely, since in the experimental regime where eV is as large as 8E T h , ac Josephson oscillations should be strongly decreased. In addition, the amplitude of the anomalies is nearly independent of the voltage value, pointing towards a fully coherent phenomenon. Figure 8 indeed shows why the quartet mechanism is immune against the decoherence brought by such voltage in a single Josephson junction. Two Andreev reflections in S 0 produce two splitted Cooper pairs, materializing in electrons and holes at energies eV ± ε in the a-branch, and A perturbative non-equilibrium Green's function [12] calculations gives the following estimate for the maximum quartet current:

eI Q ∼ -G CAR E T h sin(ϕ a0 + ϕ b0 ), (4)
where the characteristic phase dependence of the quartet mode stems from the four involved Andreev reflections, one at S a , one at S b and two at S 0 . The conductance G CAR is the Crossed Andreev conductance of a N a N S 0 N N b structure in which the electrodes S a and S b are in the normal state and at voltages ±V . By "CAR conductance" we mean the conductance formally associated to a process in which an electron from N a , at energy eV is transmitted as a hole in N b at energy -eV , and a Cooper pair is therefore split from S 0 [19]. This formula is similar to the one for a SNS junction where the pair critical current eI c ∼ G N E T h is related to the single particle normal conductance and the Thouless energy sets the coherence of Andreev reflection.

The CAR conductance can be evaluated as G CAR ∼ G N a G N b G 0 (ξs) , where G N a,b is the conductance within each normal branch of the bijunction, and G 0 (ξ s ) = (2e 2 /h)N (l e /ξ) is the normal-state conductance of a region of size ξ of the superconductor S 0 (N is the channel number). This shows that the ratio between the quartet maximum current at a bias V and the single junction critical current at zero bias is I Qmax (V )/I c (0) ∼ G CAR /G N ∼ G N /G 0 (ξ s ), which is not necessarily small. Based on measured sample parameters, we estimate this ratio to ∼ 0.1 -0.5, in fair agreement with the experiment. Notice that if eV ≪ ∆, then G CAR as well as I Q,max do not decrease with V , in agreement with the present experiment.

Conclusion

The theory and the experiment reported here open the way to a new kind of Josephson devices. Phase coherence is possible between three superconductors biased at suitable voltages, and open new channels for transport, made of correlated Cooper pairs. Such phase coherence remains to be demonstrated by interference experiments. Beyond the understanding and probing of the microscopic properties of such bijunctions (made of quantum dots or metals), future prospects include the exploration of new entanglement properties, and the production of phase-correlated microwave fields.

Figure 1 .

 1 Figure 1. A Josephson bijunction. Superconductors S a,b are biased at voltages V a,b , S 0 is grounded. The distance between the two quantum dots is comparable to the coherence length. Gates (in grey online) are figured.

Figure 2 .

 2 Figure 2. (color online)"Broad" dots regime: |ε a,b | = 6∆, Γ a,b = 4∆. Left: phase-dispersion (in blue online) of the d.c. current for the main resonances (with |n| + |m| ≤ 3), centered around the values of the phase-independent quasiparticle current (small horizontal bars). Horizontal axis is V a /V b , with V b /∆ = 0.3. Upper middle: Current I a , for the resonance 2V a + V b = 0, as a function of the phases ϕ (0) a and ϕ (0) b , showing the dependence in 2ϕ (0) a + ϕ (0) b . Lower middle: Current-phase relation I a (ϕ (0) b ) at ϕ (0) a = 0 for the quartet resonance V a + V b = 0, which shows the π-phase behavior. Right: Energy diagram for the S a D a S 0 D b S b bijunction, with a higher order diagram associated with a "sextet" current with 3 pairs emitted from S 0 , 2 into S a and 1 into S b , with V a = -V b /2.

Figure 3 .

 3 Figure 3. (color online) Currentphase relations I a (ϕ (0) b ) (red, full curve) and I b (ϕ (0) b ) (blue, dashed curve) in the quartet configuration V a = -V b , for the resonant dots regime: ε a = -ε b = 0.4∆, Γ = 0.1∆, V = 0.09 -0.775∆. Note the different y scales in the different panels.

Figure 4 .

 4 Figure 4. (color online) Upper panel: critical current I Q c a as function of the voltage in the quartet configuration, for resonant dots (same parameters as in Figure 3. Lower panel: Lowest order diagram contributing to the phasedependent MAR current.

Figure 5 .

 5 Figure 5. Experimental set-up [28]. The three macroscopic resistors have low resistance values (≃ 0.1 Ω) allowing voltage-biasing of the samples.The SEM image shows a bijunction sample with a T-shape geometry.

Figure 6 .Figure 7 .

 67 Figure 6. Differential resistance R dif f,a of a T-shape junction in the (V a , V b ) plane for various temperatures.

Figure 8 .

 8 Figure 8. Schematic diagram for Q 0 quartet production, from S 0 to S a and S b . Two Cooper pairs are split simultaneously at S 0 with one electron of each pair diffusing towards S a and S b where, under the appropriate energy condition (V a = -V b ), they recombine to create two separated Cooper pairs.

Acknowledgments

This work has been partially funded by ANR-NanoQuartet (ANR12BS1000701).

We acknowledge the Nanoscience Foundation for the PhD grant of A. H. Pfeffer and the NanoFab facility at Institut Néel-CNRS for sample fabrication.