S Samadi

H Soltanian-Zadeh

C Jutten

Cipce

Integrated Analysis of EEG and fMRI Using Sparsity of Spatial Maps

Keywords: Integration Analysis, Referenced-Based Signal Processing, Weighted Sparse Decomposition, Elastic Net, Pareto Optimization

Integration of electroencephalography (EEG) and functional magnetic resonance imaging (fMRI) is an open problem, which has motivated many researches. The most important challenge in EEG-fMRI integration is the unknown relationship between these two modalities. In this paper, we extract the same features (spatial map of neural activity) from both modality. Therefore, the proposed integration method does not need any assumption about the relationship of EEG and fMRI. We present a source localization method from scalp EEG signal using jointly fMRI analysis results as prior spatial information and source separation for providing temporal courses of sources of interest.

The performance of the proposed method is evaluated quantitatively along with multiple sparse priors (MSP) method and sparse Bayesian learning (SBL) with the fMRI results as prior information. Localization bias (LB) and source distribution index (SDI) are used to measure the performance of different localization approaches with or without a variety of fMRI-EEG mismatches on simulated realistic data. The method is also applied to experimental data of face perception of 16 subjects.

Simulation results show that the proposed method is significantly stable against the noise with low localization bias. Although the existence of an extra region

Introduction

EEG is a non-invasive technique that records scalp electrical activity generated by brain structures and it greatly helps to figure out the aspects of cognitive processes.

The EEG data is analyzed by solving its inverse problem, which looks for the localization and dynamics of cerebral activity sources. EEG reflects neural activity directly with high temporal resolution (around millisecond), but it lacks good spatial resolution. Recently, EEG-fMRI integration has been intensively addressed due to their complementary temporal and spatial resolutions. fMRI relies on the blood oxygen level dependent (BOLD) effect and measures hemodynamic signal. fMRI has high spatial resolution around millimeters. Therefore, it seems that the fMRI high spatial resolution can assist solving the EEG inverse problem.

Different EEG and fMRI integration methods have been proposed recently which can be categorized into symmetric and asymmetric groups [START_REF] Jorge | EEG-fMRI integration for the study of human brain function[END_REF][START_REF] Rosa | EEG-fMRI integration: a critical review of biophysical modeling and data analysis approaches[END_REF]. Fusion or symmetric methods use the data of both modalities as the observation variables and consider their interactions to extract the maximum information about the neural activity [START_REF] Henson | Multimodal integration: Constraining meg localization with eeg and fmri[END_REF]. These methods are complex and most of them require accurate physiological models of the interaction between neuronal activity and BOLD signals. In asymmetric methods the data from one modality is used to constrain the inference (fMRI-informed EEG [START_REF] Murta | Electrophysiological correlates of the bold signal for eeg-informed fmri[END_REF]) or as a predictor (EEG-informed fMRI [START_REF] Dong | Local multimodal serial analysis for fusing eeg-fmri: A new method to study familial cortical myoclonic tremor and epilepsy. Autonomous Mental Development[END_REF]) of the other modality.

Here, we focus on the fMRI-informed EEG group. The fMRI-informed EEG methods apply the fMRI result in the inverse problem of EEG to improve the spatial resolution of source localization. Two main challenges are related to EEGinformed fMRI method: 1) How to solve the EEG inverse problem? 2) How to apply fMRI information? Two main approaches are used for solving the inverse problem: non-parametric (imaging) and parametric (scanning) methods [START_REF] Baillet | Electromagnetic Brain Mapping[END_REF]. A brief review of these methods can be found in [START_REF] Bolstad | Space-time event sparse penalization for magneto-electroencephalography[END_REF]. Parametric or scanning methods use a small number of dipoles [START_REF] Bai | On the estimation of the number of dipole sources in EEG source localization[END_REF][START_REF] Michel | EEG source imaging[END_REF][START_REF] Shun | Efficient Dipole Parameter Estimation in EEG Systems With Near-ML Performance[END_REF], multi-poles [START_REF] Bandyopadhyay | A numerical approach for automatic detection of the multipoles responsible for ill conditioning in generalized multipole technique[END_REF][START_REF] Mosher | EEG and MEG: forward solutions for inverse methods[END_REF], or cortical patches [START_REF] Limpiti | Cortical patch basis model for spatially extended neural activity[END_REF][START_REF] Wagner | Human sleep EEG under the influence of pulsed radio frequency electromagnetic fields. Results from polysomnographies using submaximal high power flux densities[END_REF] and scan over any possible locations to find the best set of sources to represent the data. The examples of scanning methods include MUSIC [START_REF] Mosher | EEG and MEG source localization using recursively applied (RAP) MUSIC[END_REF], beamforming [START_REF] Fuchs | Anatomically constrained minimum variance beamforming applied to EEG[END_REF][START_REF] Van Veen | Minimum variance beamforming with soft response constraints[END_REF], and maximum likelihood estimation [START_REF] Inan | AR spectral analysis of EEG signals by using maximum likelihood estimation[END_REF]. Non-parametric methods are also referred to as imaging methods. In these models several dipoles with fixed locations and possibly fixed orientations are distributed in the whole brain volume or cortical surface. In imaging method, the EEG forward problem is estimated with a linear system equations. Then, using head model, estimating the sources from scalp EEG is an underdetermined and ill-posed problem, which requires prior information and constraints for converging to a unique solution.

Minimum Norm Estimation (MNE) is one of the earliest methods to solve the EEG inverse problem. MNE searches a solution with the minimum 2 norm, obtained from a closed-form equation involving a simple matrix product [START_REF] Tikhonov | Solutions of ill-posed problems. Scripta series in mathematics, Winston Tipping ME[END_REF]. This makes the estimation extremely fast. However, using the 2 norm regularization method leads to smooth solution. Therefore, it fails to recover solutions with high spatial frequencies, because the extent of active regions is often over-estimated.

Recently, sparsity has been used as an effective tool to find sparse (unique) solutions for an underdetermined system of linear equations. Sparsity means that the solution should only have a small number of non-zero elements. The 0 pseudonorm1 is the mathematical way to measure sparsity of a vector x, which is nothing but the number of non-zero components of x. In the last few years, many algorithms have been developed to solve the ill-posed underdetermined linear system when considering sparse priors [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Bruce | Block coordinate relaxation methods for nonparamatric signal denoising[END_REF][START_REF] Daubechies | Iteratively re-weighted least squares minimization: Proof of faster than linear rate for sparse recovery[END_REF][START_REF] Efron | Least Angle Regression[END_REF][START_REF] Friedman | Regularization Paths for Generalized Linear Models via Coordinate Descent[END_REF][START_REF] Gorodnitsky | Neuromagnetic source imaging with focuss: a recursive weighted minimum norm algorithm[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Tibshirani | Sparsity and Smoothness via the Fused Lasso[END_REF].

It has been proved that 1 regularization leads to a sparse solution [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and under some conditions, produces the sparsest solution [START_REF] Candès | The dantzig selector: statistical estimation when p is much larger than n[END_REF][START_REF] Donoho | For most large underdetermined systems of linear equations the minimal 1-norm solution is also the sparsest solution[END_REF][START_REF] Donoho | Compressed sensing[END_REF][START_REF] Donoho | Maximal Sparsity Representation via 1 Minimization[END_REF]. The p norm is a real-valued function which is convex when p ≥ 1. For p = 1, 1 norm is marginally convex and is the best convex estimation of the 0 norm [START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF]. Recently, some theoretical and practical studies show that the 0 regularization criterion is superior to that of 1 regularization in a sense of accuracy (Babaie-Zadeh and Jutten, 2010; [START_REF] Mancera | L0-norm-based sparse representation through alternate projections[END_REF][START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed L0 norm[END_REF]. Since 0 norm regularization leads to NP-hard problems, [START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed L0 norm[END_REF] proposed an algorithm to minimize the 0 norm by approximating 0 norm by a continuous and differentiable function. The algorithm, called smooth 0 norm (SL0), searches the optimal solution based on gradient ascent method. SL0 brings us two new opportunities. First, It is faster than the interior-point linear programming (LP) solvers, while at the same time having the same or better accuracy. Second, contrary to the previous approaches, 0 norm minimization algorithms can provide a stable estimation from noisy data, with a weaker sparsity requirement on the data than using 1 norm2 (Babaie-Zadeh and Jutten, 2010). Some studies have solved the EEG inverse problem for a sparse solution such as the self-coherence enhancement algorithm [START_REF] Yao | A Self-Coherence Enhancement Algorithm and its Application to Enhancing Three-Dimensional Source Estimation from EEGs[END_REF], the 1 norm solution (i.e., LASSO) [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], and the p norm iterative sparse source localization (LPISS) [START_REF] Xu | Marginal fisher analysis and its variants for human gait recognition and content-based image retrieval[END_REF]. Also, the combination of the 2 and 1 norms has been proposed [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF][START_REF] Gramfort | Time-frequency mixed-norm estimates: Sparse M/EEG imaging with non-stationary source activation[END_REF]. Different strategies are used to combine the 2 and 1 norms. Mixed-Norm Estimates (MxNE) algorithm minimizes the 1 norm of the 2 norm of the time series of the sources

(M xN E = M i=1 2 T j=1 |X ij | 2 ,
where X is an unknown matrix (M × T )). It also estimates the time series of the sources and their spatial maps simultaneously, but with a high computational cost. Therefore, it is used in a region of interest but cannot be applied for the whole brain. Group LASSO is used when we can partition the sparse variable into different groups [START_REF] Cassidy | Grouped l0 least squares penalised magnetoencephalography[END_REF][START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. Elastic Net uses the weighted sum of the 2 and 1 norms (λ 1 x 1 + λ 2 x 2 , where x is an unknown vector) in the regularization term to estimate the connected sparse regions.

Here, we are interested in the brain neural activities, which are sparsely localized, thus a relevant constraint is to use the sparse characteristics of the spatial distribution of the sources [START_REF] Ding | Reconstructing Cortical Current Density by Exploring Sparseness in the Transform Domain[END_REF][START_REF] Ding | Sparse cortical current density imaging in motor potentials induced by finger movement[END_REF][START_REF] Friston | Multiple sparse priors for the m/eeg inverse problem[END_REF][START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF]. We choose the 0 norm to represent sparsity and benefit its accuracy and stability against the noise using the smooth 0 norm (SL0) algorithm [START_REF] Babaie-Zadeh | On the stable recovery of the sparsest overcomplete representations in presence of noise[END_REF][START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed L0 norm[END_REF]. For using sparsity, since EEG is not potentially sparse, we have to project the EEG data to another space in which the data are represented by sparse unknown variables (see Figure 1). This space is made up of temporal sources which we are interested in, i.e. sources which are highly correlated with the paradigm of the task of interest. We call them highly correlated with interested task paradigm sources (HCTP sources). The HCTP source space is found by applying the reference-based source separation (R-SS) method proposed in [START_REF] Samadi | Reference-Based Source Separation Method For Identification of Brain Regions Involved in a Reference State From Intracerebral EEG[END_REF]. R-SS is a semi-blind source separation method which extracts the discriminating sources of two data groups, one group being related to the task of interest. By projection to the HCTP source space, we reformulate the EEG inverse problem to achieve a sparse unknown variable which is the spatial map of HCTP sources. For a given sparsity level (i.e. 0 norm), we would like to promote neural activity concentrated in neighbour pixels of a small number of cerebral areas, instead to neural activity distributed on isolated voxels.

For enforcing this sparsity, usually called block-sparsity, we will combine the 0 norm and the 2 norm like what was done in elastic net. We also use the spatial map extracted from the thresholded F-map of the fMRI data as a weight for the spatial map of the EEG source activation. It means that EEG source activity is assumed to be located in areas of significant fMRI activation with high probability.

Here, we use the same feature, extracted from EEG and fMRI, which is the spatial map of neural activity. Therefore, we did not need to use any assumption about the relationship between EEG and fMRI, which may not be generally congruent with reality.

The different spaces used in our method are illustrated in Figure 1, and defined as follows:

-Spatial map means the map of activation inside the brain.

-Temporal cerebral sources are the signals related to the cerebral sources extracted from the scalp EEG data using the source separation method.

-HCTP sources are the temporal sources (signals) highly correlated with the activity of interest: they are the outputs of the source selection algorithm.

-HCTP sources space is the space spanned by the HCTP sources, i.e. whose HCTP sources are the basis.

The main contributions of the current work are:

1. Using the same feature from both modalities which is the spatial map of the neural activity. The rest of the paper is organized as follows. The proposed method is explained in Section 2. The data protocol and the results are presented in Sections 3 and 4, respectively. Discussion and concluding remarks are brought in Sections 5 and 6, respectively.

Method

In the imaging method a grid of thousands of dipoles is used to model the forward problem of EEG as:

X = GJ + n (1)
where X is a N × T matrix that contains scalp recordings, N is the number of channels and T is the number of time points, G N ×M is the gain matrix, J M ×T contains dipole magnitudes at different time instants, n represents the noise, and M is the number of mesh vertices with, practically, N M . Under this notation, the inverse problem consists in finding an estimate of the dipole temporal courses, which are the rows of the matrix J, given the electrode positions and scalp recordings X and using the gain matrix G, which can be calculated with boundary element method (BEM).

As denoted before, the EEG inverse problem is underdetermined (due to the fewer observations than the unknown variables (since N M )) and ill-posed so, it lacks a unique solution. Therefore, some additional priors as constraints are needed to find a unique solution. The constraints are applied by adding a regularization term to the usual term of data fit, leading to the solution:

Ĵ = arg min J X -GJ + λh(J) (2)
Here, we are interested in localized brain neural activities, i.e. only a few regions of the brain that are responsible for the specific brain activity. Therefore, using the spatial sparsity as a constraint in solving the EEG inverse problem is consistent with this prior (Op de [START_REF] De Beeck | Interpreting fMRI data: Maps, modules, and dimensions[END_REF][START_REF] Simoncelli | Natural images statistics and neural reperesentation[END_REF]. The most important question is: How can we find the spatial map of neural activity to apply sparsity constraint?

We find the source spatial map in four steps. First by using a referenced based source separation (R-SS) method [START_REF] Samadi | Reference-Based Source Separation Method For Identification of Brain Regions Involved in a Reference State From Intracerebral EEG[END_REF] the time courses of the HCTP sources are estimated, secondly the forward model is projected on the HCTP source space and thirdly, the spatial map of each HCTP source is estimated using a sparse decomposition method. Finally, Pareto optimization method is used to find the active regions. From now on, as the proposed method is based on Spatial Sparsity in the Source Space, we call it "S4" method. The architecture of the proposed method is shown in Figure 2, each block of which will be explained in details in next subsections.

Source Space Estimation

Here, we use the R-SS method introduced in [START_REF] Samadi | Reference-Based Source Separation Method For Identification of Brain Regions Involved in a Reference State From Intracerebral EEG[END_REF] or rest state. R-SS uses the concept of periodic independent component analysis or common spatial pattern which can be solved using the generalized eigenvalue decomposition (GEVD) [START_REF] Koles | Spatial Patterns Underlying Population Differences In The Background EEG[END_REF][START_REF] Sameni | Multichannel Electrocardiogram Decomposition Using Periodic Component Analysis[END_REF]. Moreover, it uses the eigenvalues to estimate the number of HCTP sources (denoted i * ). The HCTP sources are the sources with the highest correlation with the reference state, which is the paradigm of the task of interest. A brief review of R-SS is explained in

SPM 1 h() Linear Model R-SS × Pareto WEN0 fMRI Data Spatial Map SM 1 h(SM ) X = GJ X s = A s S J s = B s S S =    s 1 . . . s i *    A s = G s B s B s B optimum

Projection

Fig. 2: Diagram of the proposed method (called S4). WEN0 is the regularization method, and h(SM ) = (1 -α) + αSM appendix A. The EEG recordings can be assumed as a linear combination of uncorrelated sources:

X = AS (3)
where, A is a N × N mixing matrix, a j is its j-th column, S is a N × T matrix, containing uncorrelated sources, and s i is its i-th row. Without the loss of generality the temporal sources are normalized, which means that SS = I, where I is an N × N identity matrix, and S denotes the transpose of S. We use R-SS to extract the i * HCTP sources. Therefore, Eq. ( 3) can be rewritten as:

X = [a 1 . . . a i * ] [s 1 . . . s i * ] + [a i * +1 . . . a N ] [s i * +1 . . . s N ]
We rename the HCTP part of A as A s = [a 1 . . . a i * ] and the rest is renamed as 3) can be rewritten as:

A s = [a i * +1 . . . a N ]. So Eq. (
X = X s + X s = A s [s 1 . . . s i * ] + A s [s i * +1 . . . s N ] (4)
where the dimensions of A s and A s are N × i * and N × (N -i * ), respectively.

As X and J are linearly related (if one neglects the additive noise) according to Eq. ( 1), J can also be estimated from S as:

J = BS (5)
where B is an M × N unknown mixing matrix which predicts J from S.

We can partition J in two parts estimated by HCTP sources and other sources, as done for X:

J = J s + J s = B s [s 1 . . . s i * ] + B s [s i * +1 . . . s N ] (6)
where the dimensions of B s and B s are M × i * and M × (N -i * ), respectively.

We put the HCTP sources, estimated with R-SS, in the rows of the matrix

Ŝ = [ŝ 1 • • • ŝi * ]
and call it from now on the HCTP source space.

Spatial Localization

Source Space Projection

In this step, we project the EEG forward problem to the HCTP source space by multiplying the forward model with the matrix Ŝ. Therefore, the projected version of the forward problem would be as:

X Ŝ = GJ Ŝ + n Ŝ (7)
Using Eq. ( 4) and Eq. ( 6) and the fact that sources are uncorrelated and normalized (SS ∼ I) and sources are uncorrelated with noise (n Ŝ ∼ 0), we would have:

A s I + A s0 = GB s I + GB s0 + 0 ( 8
)
where I is an i * × i * identity matrix and 0 is a zero matrix. Therefore, A s and B s are related by the following equation:

A s = GB s (9)
The i-th column of B s contains the projection of the current dipole signal of all vertices on the i-th HCTP source. Therefore, the active vertices related to the i-th HCTP source have high value in the i-th column of B s . In other words, each column of B s represents the corresponding spatial map of the HCTP sources. Eq.

(9) shows A s , the projection of the observations on the HCTP source space, is a linear combination of the spatial maps. As the brain activity is localized, most of the vertices are not active and have a small value (near zero) in B s , so B s is a sparse matrix.

Weighted Elastic Net

A newly developed variable selection method [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], called elastic net (EN), can produce a sparse model with good prediction accuracy, while encouraging a grouping effect. Experimental results and simulations have demonstrated the superiority of the elastic net over LASSO [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF].

Weighted elastic net (WEN) has also been developed in [START_REF] Hong | Weighted Elastic Net Model for Mass Spectrometry Imaging Processing[END_REF][START_REF] Jun-Tao | An Improved Elastic Net for Cancer Classication and Gene Selection[END_REF]. For a linear model

a s k = Gb s k (b s k -an M × 1 vector-is the k-th column of B s and a s k -an N × 1 vector-is the k-th column of A s ), weighted elastic net estimates are computed as: bs k = arg min b s k a s k -Gb s k + λ 2 W 2 b s k + λ 1 W 1 b s k 1
where, . and . 1 represent the 2 and 1 norms, respectively, and W 1 and W 2 are M × M weighting matrices. This optimization problem can be written as:

minimize b s k W 1 b s k 1 subject to a s k -Gb s k + λ 2 W 2 b s k < δ
Here, we replace 1 with 0 (denoted . 0 ) in the weighted elastic net and call it WEN0: minimize

b s k W 1 b s k 0 subject to a s k -Gb s k + λ 2 W 2 b s k < δ
Using the following change of variables [START_REF] Hong | Weighted Elastic Net Model for Mass Spectrometry Imaging Processing[END_REF]:

G * = (1 + λ 2 ) -1 2 G √ λ 2 W 2 a * k = a s k 0 b * k = √ 1 + λ 2 b s k
WEN0 can be written in the simplified optimization problem: The fMRI spatial map is the thresholded F-map which is normalized by its maximum value and it is registered on the mesh vertices using SPM8, which uses the value of the nearest neighbor. It is defined as SM , which is a M × 1 vector.

minimize b * k W 1 b * k 0 subject to a * k -G * b * k < δ (10)
Denoting SM i the i-th component of SM, the diagonal elements of W 1 and W 2 are initialized as:

diag(W 1 ) i = diag(W 2 ) i = 1 (1 -α) + αSM i (11)
SM provides maximum values while the EEG localization is a minimization problem (Eq. ( 10)) so the inverse of fMRI spatial map is used in the weighting matrix. In other words, a low value for a given vertex implies that it is not crucial in the sparsity of the variable and the algorithm will focus on decreasing the values of the other vertices. Equation ( 11) can be replaced by any other well-behaved function in which the weights have the inverse relationship with the fMRI spatial map. The optimum function is not known and it will be investigated in future studies.

Although, we do not know exactly the relationship between EEG and fMRI, we extract the same feature from both i.e. spatial map of the HCTP sources. b s i is the spatial map of the i-th HCTP temporal source, and it is weighted with the same feature (spatial map) from fMRI. In fact, this is very important because we are sure that we did not use any unreal assumption about the relationships between EEG and fMRI.

α can be estimated experimentally. α = 0 removes the effect of the spatial map of the fMRI data and α = 1 constrains the results strictly to the spatial map of the fMRI. Experimentally [START_REF] Liu | Spatiotemporal imaging of human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations[END_REF], the best value of the parameter α (of Eq. ( 11)) is equal to 0.9 and this value has been used by others in the fMRI-informed EEG methods. However, it is not proved theoretically [START_REF] Babiloni | Multimodal Imaging from Neuroelectromagnetic and Functional Magnetic Resonance Recordings[END_REF][START_REF] Dale | Dynamic statistical parametric mapping: combining fmri and meg for high-resolution imaging of cortical activity[END_REF]. Here, we do some experiments (see section 4.1.2) to check if the proposed value is also optimum for our experiment.

For i * times (number of HCTP sources), we solve Eq. ( 10) for b = b s k and

a = a s k , k = 1 • • • i * ,

Method Evaluation

Evaluation of the method is done in two steps. First, the method is applied to simulation data. Secondly, it is applied to experimental data (see detail in the next section).

Simulation data. Three simulation data sets are created as follows:

-The first data set includes 50 EEG-fMRI data. In this set, two sources are activated both in EEG and fMRI in the same area. This data set is produced for 6 different SNR values of the EEG data.

-The second data set includes 50 EEG-fMRI data. In this set, three sources are activated. Two of them are activated both in EEG and fMRI in the same area, and the third one, which is called fMRI extra source, is activated just in the fMRI data. In this data set, the SNR value of the EEG data is equal to 0 dB.

-The third data set includes 50 EEG-fMRI data. In this set, three sources are activated. Two of them are activated both in EEG and fMRI in the same area, and the third one, which is called EEG extra source, is activated just in the EEG data. In this data set, the SNR value of the EEG data is equal to 0 dB.

We also use this data set to show the effect of the parameter α in (11).

To compare the performance of the methods, we use two measures. The first measure is the localization bias (LB), which is defined as Euclidean distance between the mean site of the estimated activity area (μ) and the actual mean site of the simulated activity area, µ.

LB = µ -μ (12) with μ = p b p r p p b p (13)
where, r p is the MNI coordinate of a vertex and b p is its related activation level, which is the p-th row of B s (b p = (b s p ) ).

The second measure is the source distribution index (SDI), which is defined as the localization spread:

SDI = p b p r p -μ p b p (14)
The details of the simulated data are explained in section 3.1.

Experimental data. S4 has been also validated on experimental data and its group analysis is compared with the previous studies. The stability of the method against the value of α is investigated, using the idea of split-half resampling [START_REF] Strother | Evaluating fMRI preprocessing pipelines[END_REF]. In this method, the subjects are divided in two equal groups. One group is used as a train data set, while the other one is used as test data set. The results, obtained from the train group is considered as ground truth and the results of the test group is evaluated using this ground truth.

3 Data

Simulated Data

Monte Carlo modeling (sampling over randomized 50 source locations) is used to compare different methods with the proposed method. For each activation site two sets of EEG and fMRI data are created (details are given in the two next subsections). 50 random vertices are selected as the center of their active regions.

The sources are assumed to be spread, so a smoothing Gaussian filter is used to model the spread sources.

Here, we simulate the EEG and fMRI data after the preprocessing step. It means that for example in fMRI, slice timing is not simulated and in EEG, the mean of the trials after preprocessing is simulated.

Sources are located in the grey matter. Grey matter is extracted from the tessellation of the real structural MRI of the first subject of our experimental data set. To be more realistic, we performed the source localization on the experimental EEG data and the spatial map is considered as an ROI (region of interest) in the simulated data set. We have put the sources in the ROI.

The experiment has two states: task and control. Onsets of task and control states are set as the onsets of the experimental data for face and scrambled states, respectively. All the parameters of imaging methods are chosen as the experimental data, like the electrode positions and sampling rate in EEG and TR and image dimension in fMRI.

fMRI

We use the package neuRosim in R (R Core Team, 2012) for fMRI simulation, because the data generation in neuRosim is accurate and fairly fast and the parameters can be adjusted by the user [START_REF] Welvaert | neu-Rosim: An R Package for Generating fMRI Data[END_REF]. A double-gamma function which is a linear combination of two gamma functions is used to model HRF [START_REF] Friston | Eventrelated fMRI: characterizing differential responses[END_REF][START_REF] Glover | Deconvolution of impulse response in event-related BOLD fMRI[END_REF]. Fifty fMRI simulations are created with noisy measurements with a fixed SNR= 5dB. In fact, the main known effect of low SNR in the result of our method is due to the number of false positive or false negative regions. Therefore, we evaluate the efficacy of our method in the presence of fMRI invisible or extra sources, at a fixed SNR, defined as:

SN R = 20 log( S σ N ) ( 15
)
where S represents the average magnitude of the signal, and σ N stands for the standard deviation of the noise [START_REF] Krüger | Physiological Noise in Oxygenation-Sensitive Magnetic Resonance Imaging[END_REF]. The noise is a mixture of different random signals: 1) Rician system noise; 2) temporal noise of order 1;

3) low-frequency drift; 4) physiological noise; and 5) task-related noise [START_REF] Krüger | Physiological Noise in Oxygenation-Sensitive Magnetic Resonance Imaging[END_REF]). The background is set to the mean image of the experimental data.

In each fMRI data two sources are estimated: one is related to the task state and the other to the control state (see Figure 3). The center of the spatial sources related to the task state is chosen randomly and the sources are extended with a radius of 10 mm. The control source is located at (-10,-29,55) mm in the Montreal Neurological Institute (MNI) coordinates and is extended as a sphere with a radius of 7 mm and a fading ratio of 0.01. The task sources are located at least 6 mm far from the control source.

Extra source in fMRI. For evaluating the effect of extra fMRI sources which cannot be detected in the EEG data, another data set is produced by appending another active region in the fMRI data but not in the EEG data. The extra region is supposed to be located in a fixed area and the other source which is active both in EEG and fMRI in the same area (common source) is chosen randomly in the ROI at least 8mm far from the extra source in fMRI.

EEG

For each fMRI data, the related EEG data is generated using the linear forward model (Eq. ( 1)). The gain matrix is calculated using the structural MRI of a real subject and the position of the electrodes in experimental data acquisition.

The time courses of the sources are generated as an additive source problem in which two main sources are distinguished, namely (1) the activation caused by an experimental design, and (2) the background. The experiment is defined by two states, task and control trials. The details are as follows.

-Activation site. We assume 50 random sites to put a dipole related to the activation of interest (task). These random sites are the same for EEG and fMRI. We also assume a fixed site for a dipole related to the control state (like in fMRI it is located at (-10,-29,55) mm in the MNI coordinates).

-Activation time courses. For simulating the activation time courses, we extract EEG source time courses from the experimental data. To extract the task and control time courses, the R-SS method is used on the experimental data. To do this, alternatively, the task interval (faces) or control state (scrambled faces) is assumed to be the reference state, for estimating the task temporal course or the control temporal courses, respectively.

-Background time course. It is extracted from experimental data. It can be found by extracting the common sources between the task and control trials, X 1 and X 2 , respectively. Common source extraction is done by canonical correlation analysis (CCA) [START_REF] Hotelling | Relations between two sets of variates[END_REF]. We use the first common source as an estimation of the background EEG signal. By back reconstruction, we compute the EEG background in the electrode space.

-Time courses of the electrodes. To compute the electrode time courses, we use the forward linear model (Eq. ( 1)). To estimate the gain matrix, we use the structural MRI of a real subject with its fiducials and electrode positions. For this purpose, we use BEM with 8196 nodes developed in SPM8

(http://www.fil.ion.ucl.ac.uk/spm). As EEG is sensitive to the signal of the Pyramidal neurons of the cortex, which are perpendicular to the cortical surface, we assume that the dipoles are on the cortical grey matter, where each dipole is perpendicular to the cortical sheet. We put the task and control time courses in the related rows of the matrix J. We put the estimated activation sources in the predefined sites and use a Gaussian filter with full width at half maximum (FWHM) of 8 mm to estimate the activation of the neighboring mesh vertices. After computing the effects of activations on the electrode space, we add the background signal to the electrodes with different SNR values in the range of [-20, 5] dB with a step size of 5 dB.

Extra source in EEG. For evaluating the effect of extra EEG sources which cannot be detected in the fMRI data, another active region in the EEG data is added with a weak signal power which is not detected in the fMRI data.

Experimental Data

Experimental data is taken from ftp://ftp.mrc-cbu.cam.ac.uk/personal/rik.

henson/wakemandg_hensonrn/. A brief summary of the acquisition characteristics, and preprocessing steps are brought in this section. More details can be found in [START_REF] Henson | MEG and EEG data fusion: Simultaneous localisation of face-evoked responses[END_REF][START_REF] Henson | A Parametric Empirical Bayesian framework for the EEG/MEG inverse problem: generative models for multisubject and multimodal integration[END_REF].

Participants

The participants are sixteen healthy young adults (eight female). The study protocol is approved by the local ethics review board (CPREC reference 2005.08).

Each subject performs two days experiments; once for concurrent MEG + EEG and once for fMRI + sMRI. The protocol is the same each day, but the fMRI and EEG data are not recorded simultaneously.

Stimuli and Tasks

The paradigm is similar to that used previously under EEG, MEG, and fMRI [START_REF] Henson | MEG and EEG data fusion: Simultaneous localisation of face-evoked responses[END_REF][START_REF] Henson | A Parametric Empirical Bayesian framework for the EEG/MEG inverse problem: generative models for multisubject and multimodal integration[END_REF]. A central fixation cross (presented for a random duration of 400-600 ms) is followed by a face or scrambled face (presented for a random duration of 800-1000 ms), followed by a central circle for 1700 ms.

EEG Acquisition

The EEG data is recorded simultaneously with MEG with a VectorView system (Elekta Neuromag, Helsinki, Finland). EEG data is recorded with 70 Ag-AgCl electrodes with the sampling rate of 1.1 kHz. The low-pass filter with the cut-off frequency of 350 Hz is applied on EEG data.

An elastic cap (EASYCAP GmbH, Herrsching-Breitbrunn, Germany) is used according to the extended 10-10 and a nose electrode is used as the recording reference. Vertical and horizontal EOG (and ECG) are also recorded.

fMRI + sMRI Acquisition

Structural MRI (sMRI) is a T1-weighted MPRAGE-MRI scan, which was acquired for each subject with voxel size 1 × 1 × 1 mm, on a 3T Trio (Siemens, Erlangen, Germany) with TR= 2250 ms, TE= 2.99 ms, flip-angle= 9 • and acceleration factor= 2.

The fMRI volumes comprise 33 T2-weighted transverse echoplanar images (EPI) (64 × 64, 3 mm × 3 mm pixels, TE = 30 ms) per volume, with blood oxygenation level dependent (BOLD) contrast. EPIs comprise 3 mm thick axial slices taken every 3.75 mm (3 mm thick with a 0.75 mm gap), acquire sequentially in a descending direction. 195 volumes were collected continuously with a repetition time (TR) of 2.25s. The first six volumes were discarded to allow for equilibration effects.

EEG Preprocessing

The steps of preprocessing are as follows:

-Epoching from -500 to +1000 ms poststimulus onset.

-Down-sample to 250 Hz (using an antialiasing low-pass filter with a cutoff frequency of approximately 100 Hz).

-Reject the epochs with the EOG exceeded 150 µV (number of rejected epochs ranged from 0 to 47 across participants, median=5).

-Re-reference the EEG data to the average over non-bad channels. The number of bad channels, which are rejected, varies between 0 to 7, according to the participants.

sMRI Processing

Structural MRI is used to model the forward problem. sMRI images of each participant are segmented and spatially normalized to an MNI T1 template brain in Talairach space. The inverse of the normalization transformation is then used to warp a cortical mesh from a template brain in MNI space to each participant's MRI space. The mesh has 8196 vertices (4098 per hemisphere) with a mean intervertex spacing of 5 mm. The lead-field matrix is then calculated for the dipoles at all points in the cortical mesh, oriented normal to the mesh. All is done in SPM.

fMRI Processing

The processing steps are performed with SPM8 (http://www.fil.ion.ucl.ac.uk/ spm). General linear model (GLM) is used to model BOLD response and maximum likelihood estimation is used to estimate the parameters. A statistical parametric map (SM) of the F-statistic that compare faces against scrambled faces was thresholded for p < 0.05 (family-wise error corrected across the whole-brain) and the regions of at least 10 contiguous voxels are considered as active regions.

This analysis produces clusters in the fusiform and occipital face areas (FFA and OFA, respectively).

Results

The efficiency of S4 is compared with two other methods. The first method is the multiple sparse priors (MSP) [START_REF] Friston | Multiple sparse priors for the m/eeg inverse problem[END_REF] with or without the fMRI spatial map as prior information. MSP is implemented in SPM8 (http://www.fil.

ion.ucl.ac.uk/spm). We also compare S4 with sparse Bayesian learning (SBL)

with or without the fMRI spatial map as prior information [START_REF] Wipf | A unified Bayesian framework for MEG/EEG source imaging[END_REF]. Among the integration methods MSP [START_REF] Friston | Multiple sparse priors for the m/eeg inverse problem[END_REF] and SBL (Tipping, 2001) are the most similar methods to the proposed one, as they use the spatial sparsity of the activated regions. In these methods the active sources can be correlated and any extra information which gives information about the covariance matrix of the sources can be used e.g. the fMRI results. In these methods the information about the experiment is not used, while we have used it in the proposed method.

Simulated Data

Spatial Accuracy

To compare the spatial accuracy of different methods, mean values and standard deviations of LB (Eq. 12) are calculated for all methods (see Figure 4). Also, mean values and standard deviations of SDI (Eq. 14) are represented in Figure 5. It is the result of using fMRI information directly to solve the EEG inverse problem and using the 0 norm. In addition, S4 is more robust to noise. This is because the reference-based source separation method removes the common sources between the two states and thus, it removes the background effect as much as possible.

Figure 5 shows the mean source distribution index (SDI) of different methods in different SNR values. The mean SDI of the actual source is 15mm. Therefore, the results of S4 and SBL are more accurate. In contrary, MSP estimates sources with a too small SDI.

Effect of parameter α

Extra Source in fMRI. For SN R EEG = 0 dB, repeating the Monte Carlo simulation, using the second simulated data set, shows that the results would not change even for α = 0.9 when there exists an extra source in fMRI, which is not detectable in EEG. In this case the mean of LB is equal to 18 mm, and no extra region has been found. It means that our method will miss the regions which are in fMRI 11)).

To see the effect of α, we have repeated the experiment for α = 0.1, • • • α = 0.9.

Figure 6 shows LB for three different values of α calculated for the common source.

Table 1 shows the percentage of the EEG extra source detection in the 50 EEG-fMRI data sets, and Table 2 shows the percentage of the common source detection in the 50 EEG-fMRI data sets. Detection happens when S4 finds a region less than 20 mm far from the EEG extra source.

Although α = 0.9 is the optimum value, which is suggested in [START_REF] Liu | Spatiotemporal imaging of human brain activity using functional MRI constrained magnetoencephalography data: Monte Carlo simulations[END_REF],

in our experiments α = 0.7 is the optimum value. Tables 1 and2 show that for α = 0.7, we have the highest detection percentage of the extra source while the common source is estimated correctly in all experiments. We have applied MSP and SBL on the same data set, in which there exists an extra EEG source. These methods did not converge to the real solution and they fail to estimate even the common source site.

Experimental Data

Spatial Sources

The maps of spatial sources for group analysis of 16 subjects, using MSP and S4, are shown in Figures 7 and8, respectively. A zoom of the spatial map extracted by S4 algorithm is shown in Figure 9. For group analysis, the EEG data were analyzed for each subject and SPM is used to do group analysis. One-sample ttest is used and the t-statistics is thresholded for p < 0.05 and the regions of at least 10 contiguous voxels are considered as active regions.

The other methods, MSP and SBL, with fMRI priors could not find a common region between different subjects as their solutions are too sparse, and they did not consider any grouping constraint.

Effect of parameter α

To evaluate the stability of S4, we count the false positive and false negative regions using the idea of split-half resampling [START_REF] Strother | Evaluating fMRI preprocessing pipelines[END_REF]. For this purpose, we partition the EEG data into two non-overlapping groups each with 8 subjects.

For each group, the spatial map is extracted using S4. One group is assumed as the ground truth and the other as the test one. in the test data is at least 20 mm far from all regions in the ground truth, it would be false positive. In reverse, if any region in the ground truth is at least 20 mm far from all regions in the test data, it would be false negative. 20 mm is chosen because the best resolution of MSP in the simulation was around 20 mm (see Figure 4) and it is assumed that the regions that are less than 20 mm far from each other are acceptable as close regions. This test is repeated 20 times for 20 different partitioning. The results for two different values of α (0.9 and 0.7) are shown in table 3. For α = 0.9, no false positive region is detected in all experiments and only two false negative regions occurred in all experiments. For α = 0.7, two false positive regions are detected in all experiments and eight false negative regions occurred in all experiments. The number of regions detected in each experiment is between 1 and 5 (mean=3). The results show that α = 0.9 is more stable against subject changes which was expected, because we have used the fMRI group analysis as the constraint, and the results of the EEG group analysis will be more restricted to the fMRI results, with higher value of α. Therefore the search space in different subjects is limited to the same space (fMRI group analysis spatial map) and the variability between subjects would be ignored in higher value of α. Table 3: Effect of α in group analysis.

α for the train group 0.7 0.9 α for the test group 0.7 0.9 number of false positive regions 2 0 number of false negative regions 8 2

Discussion

Although EEG-fMRI integration is intensively investigated in recent years, some challenges still exist as open problems. The most important challenge is the unknown relationship between EEG and fMRI data. In this study, we investigate the EEG inverse problem from a new point of view. By projecting the EEG inverse problem to HCTP source space, we accurately associate the temporal space to the spatial space.

The projection basis used in S4 (HCTP sources) directly results in the spatial maps of interest, which are assumed to be sparse. Therefore, we can use sparse decomposition methods: in this work, the 0 norm regularization leads to high accuracy and stability against the noise.

After projection, the unknown variables are HCTP spatial maps, so we could use the spatial map of activation, extracted from fMRI, in the weighting matrix.

Therefore, we solve the integration inverse problem using the same feature ex-tracted from fMRI and EEG. This strategy was a remedy for overcoming (without assumption or estimation) the unknown relationships between EEG and fMRI.

However, if we have a more accurate information about the relationship between EEG and fMRI data, it may improve the results. Especially, we could optimize the weighting matrices, instead of choosing an arbitrary function to use the fMRI spatial map in the weighting matrix.

Another challenge in EEG-fMRI integration is the mismatch between EEG and fMRI sources. The proposed method has a good stability against an extra region in fMRI. Although it was not capable of detecting EEG extra regions, it achieves an acceptable detection rate of common sources between EEG and fMRI in the presence of mismatches.

Without definitely solving the open questions of EEG-fMRI integration, experiments show the S4 method achieves a step forward to use fMRI spatial maps in the EEG inverse problem with more feasibility.

The results are discussed in details as follows.

Simulated Data

The results of the method on the simulated data and its comparison with two methods, MSP and SBL, show that the spatial maps estimated by S4 and SBL, which minimize the 0 norm, have smaller mean and variance of localization bias (LB) for all SNR values, compared to results obtained with MSP (Figure 4). SBL and S4 also estimate the mean of source distribution index (SDI) more accurately, while S4 has less variance of SDI (Figure 4). Therefore, SBL and S4 are more accurate and stable against the noise.

Invisible sources either in fMRI or EEG result in an inconsistent constraint (here, it is a weighting matrix that contains fMRI information) in the EEG inverse problem. In this case, MSP and SBL fails even to estimate the common sources between EEG and fMRI, while S4 correctly estimates the common sources.

The activation regions which are missing in fMRI data may be related to neural activities which have a short activation time or involve a few numbers of neurons that cannot produce a detectable increase of cerebral blood flow [START_REF] Babiloni | Multimodal integration of EEG, MEG and fMRI data for the solution of the neuroimage puzzle[END_REF]. Missed regions can also be due to the analysis method. The result shows that the effect of missed regions in fMRI is dependent on the balancing variable α in Eq. ( 11). With α = 0.9 fMRI invisible sources are also missed in the EEG-fMRI result, but with α = 0.7, the EEG extra source could be detected in 54 percent of our experiments. Figure 6 and Tables 1 and2 show that, by decreasing the α value, the probability of the extra source detection would increase while the mean of LB (localization bias) would not change much. The α parameter can be considered as a user-defined parameter for the method. User can decrease the α value to see if there exists any extra region in EEG. However, the optimum value of α is an open question by now.

Some brain cells (stellate cells) and cortical regions (thalamus) produce closed electromagnetic field. This field cannot be detected on the scalp electrodes, although their activation requires high blood flow. This would also happen when "sources of activation cancel each other according to the orientation disparity of the cortex" [START_REF] Irimia | Source cancellation profiles of electroencephalography and magnetoencephalography[END_REF]. Therefore, these neuronal populations present high metabolism requirements that can be detected by the fMRI technique, while at the same time they are invisible in the EEG modality [START_REF] Babiloni | Multimodal Imaging from Neuroelectromagnetic and Functional Magnetic Resonance Recordings[END_REF]. We have evaluated the algorithm S4 when it exists an extra source in fMRI that is invisible in EEG. In this case, S4 increases the mean of the localization bias from 10 mm to 18 mm, but again, it outperforms the other methods, which fail to estimate the common source in the existence of an extra source.

Experimental Data

S4 algorithm is applied on experimental EEG and fMRI data of a face perception experiment. Spatial sources are shown in Figure 8. These activation regions are located in the fusiform face area (FFA) and the occipital face area (OFA). These clusters are in general agreement with previous studies reporting a similar contrast of faces versus scrambled faces [START_REF] Henson | MEG and EEG data fusion: Simultaneous localisation of face-evoked responses[END_REF][START_REF] Henson | A Parametric Empirical Bayesian framework for the EEG/MEG inverse problem: generative models for multisubject and multimodal integration[END_REF]. Previous studies investigated the role of these two regions in face perception and they proposed that the OFA is the first stage in a hierarchical face perception network in which the OFA represents facial components (e.g., eyes, nose, and mouth) [START_REF] Pitcher | The role of the occipital face area in the cortical face perception network[END_REF], while the FFA encodes both the facial components and the configural information of faces (i.e., spatial relation among features) [START_REF] Zhang | The fusiform face area is engaged in holistic, not parts-based, representation of faces[END_REF]. The most important factor here is that this activation area should be similar in all participants.

In other words, the activation site should have little spatial variation in different participants. The results show that S4 algorithm can find a common region between subjects, while MSP and SBL are not able to estimate a common region between different subjects because their solutions are sparse but not block-sparse.

In fact, we have counted the false positive regions and false negative regions using the split-half resampling algorithm proposed in [START_REF] Strother | Evaluating fMRI preprocessing pipelines[END_REF]. The results show high stability or low spatial variation over different participant groups with α = 0.9 in which the results are strictly biased with fMRI spatial map. For α = 0.7, more false positives and false negatives are detected, but again it is acceptable.

Conclusion

We propose a new method, called S4, to integrate scalp EEG recording and fMRI data for the identification of the active regions of the brain with high spatial resolution. Using R-SS, a semi-blind source separation method, we estimate the spatial filter that maximizes the power of reference over non-reference states and provides temporal sources with high correlation with the paradigm of the task of interest (called HCTP sources). The EEG signal was projected to the HCTP source space. The inverse problem of EEG in the HCTP source space can be considered as a sparse decomposition problem. The unknown variables are the spatial maps of the HCTP sources. As we focus on the localized sources, their spatial maps are sparse. For promoting block-sparsity, we extend the elastic net idea by combining 0 and 2 norms. This algorithm, called WEN0, takes into account the fMRI information through weighting matrices related to the F-map of fMRI extracted with GLM.

In more details advantages of the proposed method are as follows.

-The simulation results show that the localization quality of S4 and SBL, measure with localization bias (LB) and source distribution index (SDI), is 50 percent better than MSP. The results of S4 is significantly stable against the noise.

-The stability of S4 against false positive and false negative regions in the prior information is investigated. S4 can estimate the common source between EEG and fMRI in the existence of an extra region and missed region in the fMRI data, while both SBL and MSP fail to estimate the sources.

-Using Pareto optimization helps us in avoiding additional statistical methods and thresholding.

The drawback of S4 is that the final result is dependent on the first source separation step and any mistake in the first step will propagate to the rest. However, we did not obtain any false discovery in the simulation experiments, even for the lowest SNR value.

In future works, S4 may be applied to simultaneous EEG and fMRI data.

Defining different weighting matrices for the 2 and 0 norms may be investigated and also, the optimized value of α, which controls the contribution of fMRI spatial map in the weighting matrix, can be defined using the data.

ranking of the estimated temporal sources, according to their resemblance to the reference activation state. Here, we propose a simple method to define the number of sources. Remember that the simplicity of this method is due to two facts; 1) the usage of experiment information in the source separation step, and 2) a special interpretation of the eigenvalues in the source selection step, which will be explained below.

After obtaining the discriminative sources (s i ) between the reference and nonreference states, sources are ranked according to their similarity to the reference state. Now, we need to select the sources, which are similar enough to the reference, for being considered as belonging to the reference class. To this end, we propose the following procedure.

The probability of the reference class (ω 1 ) membership is calculated as follows:

p(s i ∈ ω 1 ) = λ i max(λ j=1,...,N )

where λ i , i = 1, . . . , N indicate the eigenvalues or the diagonal elements of Λ in (18). For the classification of the sources, the error probability using Bayes rule is defined as:

p error = N j=1
p(s j ∈ ω 2 |ω 1 )p(ω 1 ) + N j=1 p(s j ∈ ω 1 |ω 2 )p(ω 2 ) where p(s j ∈ ω 2 |ω 1 ) = 1 -p(s j ∈ ω 1 ) and p(s j ∈ ω 1 |ω 2 ) = 1 -p(s j ∈ ω 2 )

as ω 1 and ω 2 constitute a partition. p(ω 1 ) and p(ω 2 ) are the prior probabilities of the reference and non-reference classes, respectively. We remind that GEVD sorts the separated sources in the decreasing order of similarity with respect to the reference. Therefore, if we assume that only the first i sources belong to the reference class (and consequently the N -i others belong to the non-reference class), then the total error probability (false positive plus false negative errors) can be written as follows:

p error (i) = i j=1 p(s j ∈ ω 2 |ω 1 )p(ω 1 ) + N j=i+1 p(s j ∈ ω 1 |ω 2 )p(ω 2 ) (20)

where p(ω 1 ) = i N and p(ω 2 ) = N -i N . Thus, the minimum of p error (i) provides the number i * of the sources in the reference class, i.e., i * = arg min i p error (i).

Appendix B Pareto Optimization B s is a M × i * matrix whose columns present the contributions of the related sources in each mesh vertex. All the i * sources are involved in the neural activation.

Therefore, all of them are important, and if we use the sum or the square sum of the contributions the weak sources would be hidden in the shadow of the sources with higher power. The remedy of this problem is to use multi-objective optimization method, called Pareto method [START_REF] Deb | Multi-Objective Evolutionary Algorithms: Introducing Bias Among Pareto-Optimal Solutions[END_REF]. Therefore, we are able to take into account the individual effect of each source. Pareto optimization finds the optimum without any threshold, that obviates the need for arbitrary user-defined threshold. 3. If i is not equal to M go to step 2.

Fig. 1 :

 1 Fig. 1: The different spaces and their relationships.

  2. By projecting the forward problem of EEG onto the HCTP source space, we can reach enough sparsity to use sparse decomposition methods. 3. Promoting solutions characterized by the concentration of active voxels in a few regions (block-sparsity) using the 0 norm combined with the 2 norm. 4. Using fMRI spatial map as a weighting matrix in the penalty term to solve the inverse problem of EEG. 5. Exploiting Pareto optimization to find the final results.

  Using the fast algorithm proposed by[START_REF] Babaie-Zadeh | Weighted Sparse Signal Decomposition[END_REF][START_REF] Mohimani | A fast approach for overcomplete sparse decomposition based on smoothed L0 norm[END_REF] to solve Eq. (10), we estimate b s 1 , • • • , b s i * which are the columns of B s . Stability against the noise and convergence of the algorithm are discussed in details in (Babaie-Zadeh and Jutten, 2010). In our problem, a s k (an N × 1 vector) is the k-th column of A s which is the projection of the EEG recordings onto the k-th HCTP source, and the weighting matrices W 2 and W 1 are N × N diagonal matrices initialized by the inverse of the spatial map extracted from the fMRI data.

Fig. 3 :

 3 Fig. 3: Sample activation sites in simulated data (the figure is generated with xjView toolbox (http://www.alivelearn.net/xjview)).

Fig. 4 :

 4 Fig. 4: Mean of localization bias over 50 random activation sites compared for five methods (MSP, MSP+fMRI, SBL, SBL+fMRI, S4).

Figure 4

 4 Figure 4 shows the mean localization bias of different methods in different SNR values. Comparing the results of MSP and SBL when fMRI information is used or not, shows that using fMRI information results in better resolution especially in lower SNR values. The figure also shows that SBL and S4 have the same accuracy.

Fig. 5 :

 5 Fig. 5: Source distribution index over 50 random activation sites compared for five methods (MSP, MSP+fMRI, SBL, SBL+fMRI, S4).

Fig. 6 :

 6 Fig. 6: Localization bias in the presence of an extra EEG source as a function of α.

  Fig. 7: Activation sites in the experimental data of 16 subjects by MSP (the figure is generated with xjView toolbox, http://www.alivelearn.net/xjview, and scale unit=1).

Fig. 8 :

 8 Fig. 8: Activation sites in the experimental data of 16 subjects by the proposed method (the figure is generated with xjView toolbox, http://www.alivelearn.net/xjview, and scale unit=1).

Fig. 9 :

 9 Fig. 9: Zoom of the activation sites in the experimental data of 16 subjects by S4 (the figure is generated with xjView toolbox, http://www.alivelearn.net/xjview).

A

  multi-objective optimization problem, in Pareto sense, has the following form: maximize (b s i ) for i = 1, • • • , M subject to (b s i ) ∈ P ⊂ i * (21) It then consists of i * objective functions that are aimed to be maximized simultaneously. (b s i ) is the i-th row of the matrix B s which shows the contributions of the sources at the ith vertex. From the geometrical point of view, each (b s i ) can be considered as a point in a i * -dimensional space. In Pareto optimization the non-dominated points should be chosen as the optimum points (Deb, 1999): a point is non-dominated if either it dominates the others, or there is no other point dominating it. Point (b s i ) dominates point (b s k ) , if ∀l, b s i (l) ≥ b s k (l), and ∃l * , b s i (l) > b s k (l * ) 3 , where b s i (l) is the l-th component of the vector b s i . The set of all non-dominated points is called non-dominated layer. Let us consider M i * -dimensional decision vectors, (b s i ) , as M points in the search space P . The non-dominated layer, denoted by D(P ), is obtained using the following Pareto optimization algorithm (Deb, 1999): 1. Initialize D(P ) with the first point (i = 1) with the value of (b 1 ) . This can be any point. 2. Choose a new point (i = i + 1): (a) If any node in D(P ) dominates point i go to step 3. (b) Else add point i to D(P ) and remove any points of D(P ) that point i dominates.

  which are the columns of B s and A s , respectively. Therefore, we estimate i * spatial maps related to the i * HCTP sources derived from R-SS. Now, we should investigate the HCTP spatial maps and extract a unique spatial map. To this goal, for finding the activation site or the place of optimum coefficients in B s , we use Pareto optimization (see appendix B) .

Table 1 :

 1 Detection percentage of the extra EEG source.

	α	0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
	Detection (%)	58	58	52	56	52	46	54	20	10

Table 2 :

 2 Detection percentage of the common EEG source.

	α	0.1 0.2 0.3 0.4 0.5	0.6	0.7	0.8	0.9
	Detection (%)	36	46	46	54	68	100 100 100 100

The 0 pseudo-norm does not satisfy the mathematical definition of a norm, However, in the following, we simply say 0 norm.

i.e. with data less sparse than required with algorithm based on 1 norm.

Acknowledgement

This work has been partially funded by the project CHESS, 2012-ERC-AdG-320684. The authors gratefully acknowledge Daniel Wakeman and Richard Henson for providing the actual data sets.

Appendix A Referenced-based Source Separation

The objective is to extract the sources related to a reference activation model.

We consider two states, denoted C 1 and C 2 , which correspond to the reference and non-reference activations, respectively. Denoting T , = 1, 2, the set of time samples related to each state, we can build the corresponding segment matrix, X ∈ R N ×M .The correlation matrix of data for each state can be estimated as:

The spatial filters, W (whose columns are generically denoted w), for which the temporal sources, S = W X have maximum similarity with the reference activation state, i.e., maximum variance in the reference state compared to the other state, is computed as:

Solving ( 17) leads to generalized eigenvalue decomposition (GEVD) of ( R

Using W, the spatial patterns, A = (W ) -1 , and the temporal sources, S = W X, are extracted. The maximum eigenvalue in ( 18) is related to the maximum power ratio in (17). We rank the eigenvalues in decreasing order. This implies


Fabrice Kordon
email: fabrice.kordon@lip6.fr

Hubert Garavel

Lom Messan Hillah

Emmanuel Paviot-Adet

Loïg Jezequel

César Rodríguez

Francis Hulin-Hubard

MCC'2015 -The Fifth Model Checking Contest

The Model Checking Contest (MCC) is an annual competition between software tools that verify concurrent systems using state-space exploration techniques, either explicit-state or symbolic. The present article provides a comprehensive account of the 2015 edition of the MCC. The principles of the contest are described, together with its underlying software infrastructure. The tools that competed in 2015 are listed and the results of the contest are summarized.

Goals and Scope of the Model Checking Contest

For more than a decade, one has seen the emergence of software contests that assess the capabilities of verification tools on complex benchmarks, so as to identify which theoretical approaches are the most fruitful ones in practice, when applied to realistic examples. Notable events that have significant impact on the involved communities include: the SAT competition (nine editions since 2002), the Satisfiability Modulo Theories Competition (ten editions since 2005), the Hardware Model Checking Contest (eight editions since 2007), the Verified Software Competition (four editions since 2010), the Rigorous Examination of Reactive Systems Challenge (five editions since 2010), the Timing Analysis Contest (one edition in 2011), and the Competition on Software Verification (four editions since 2012). The existence of long-lasting events is a clear indication of interest and usefulness.

The Model Checking Contest (MCC for short -five editions since 2011) belongs to this family of scientific contests. It aims at evaluating model-checking tools that analyze formal descriptions of concurrent systems, i.e., systems in which several processes run simultaneously, communicating and synchronizing together. Examples of such systems include hardware, software, communication protocols, and biological models. So far, all editions of the MCC used Petri nets to describe such systems, but there is a longterm goal to open the contest to other model-checking tools not primarily based on Petri nets.

Since the first edition MCC'2011 within the context of a workshop associated to the Petri Nets conference, the contest team has attracted key people with diverse knowledge, who are actively contributing to the selection of benchmark models, the auto-mated generation of temporal-logic formulas, and all technical aspects related to performance measurements and tool assessment. Today, the MCC team gathers scientists from CNRS, Inria, Université Pierre & Marie Curie, Université Paris Ouest Nanterre La Défense, Université Paris 13, and Université de Nantes.

The present paper reports about the fifth edition MCC'2015, which was organized in Brussels as a satellite event of the 36 th International Conference on Application and Theory of Petri Nets and Concurrency. This is the first publication devoted to the MCC since its origins, and it will be followed by detailed presentations of the tools that reached the top-three podiums in June 2015. The paper is organized as follows. Section 2 presents the collection of benchmarks on which the participating tools are assessed. Section 3 explains how temporal-logic formulas are automatically generated for the contest examinations. Section 4 details how tools are executed and how information is retrieved from their execution. Section 5 lists the participating tools, and Section 6 summarizes the results of the contest. Finally, Section 7 gives some concluding remarks and suggests desirable enhancements for future MCC editions.

Collected Benchmarks

All tools participating in a given edition of the MCC are evaluated on the same benchmark suite, which is updated every year. The yearly edition of the MCC starts with a call for models inviting the scientific community at large (i.e., beyond the developers of the participating tools) to propose novel benchmarks that will be used for the MCC. The benchmarks obtained this way are merged with those of the former years to form a growing collection (continuously expanded since 2011) that gathers systems from diverse academic and industrial fields: software, hardware, networking, biology, etc. This collection of benchmarks1 is a perennial result of the MCC organization work -should the contest halt, the collection would remain available to the scientific community. The usefulness of this collection is already witnessed by over thirty scientific publications 2 .

Models and Instances. The collection of MCC benchmarks consists of models, each corresponding to a particular academic or industrial problem, e.g., a distributed algorithm, a hardware protocol in a circuit, a biological process, etc. Models may be parameterized by one or a few variables representing quantities such as the number of agents in a concurrent system, the number of message exchanged and the like. To each parameterized model are associated as many instances (typically, between 2 and 30) as there are different combinations of values considered for the parameters of this model; each non-parameterized model has a single associated instance. Each instance corresponds to a Petri net that will actually be submitted to all participating tools. Models come from diverse sources: 34% of the models originate from high-level colored Petri nets, which are then automatically unfolded to P/T nets; 14% of the models have been produced from formal descriptions written in LOTOS [START_REF] Iso/Iec | LOTOS -A Formal Description Technique Based on the Temporal Ordering of Observational Behaviour[END_REF] or more recent process-calculi languages that translate to LOTOS, from which Petri nets can be 1. Accumulation of models and instances over the years (supplemented with the 2016 data) automatically generated using the CADP toolbox [START_REF] Garavel | CADP 2011: A Toolbox for the Construction and Analysis of Distributed Processes[END_REF]; the remaining 52% have been manually specified as genuine P/T nets, e.g., using a graphical editor. Table 1 illustrates the growth of the collection since the first edition of the MCC; the notion of NUPN (Nested-Unit Petri Net) is discussed below.

All instances in the MCC benchmark collection are provided as PNML files [START_REF]High-level Petri Nets -Part 2: Transfer Format. International Standard[END_REF]. After each annual call for models, the MCC model team examines all the files received to check their conformance to the PNML standard and, if needed, normalize their contents to avoid some known ambiguities of PNML -for instance, to give the same value to the id and name attributes of each place and each transition (thus ensuring unique naming conventions across all the participating tools), and to replace multiple occurrences of arcs between the same place and the same transition by a unique arc with multiplicity greater than one.

Forms and Properties. Each MCC model is described in a form, which is a two-or three-page PDF document giving a high-level description of the model: origin, functionality, bibliographic references, and graphical representation (if available). The form defines the model parameters, if any, and their corresponding values, as well as the size (number of places, transitions, and arcs) of each instance of the model. The form also lists essential properties of the model, both structural (e.g., strongly connected, loop free, etc.) and behavioural ones (e.g., safeness, liveness, reversibility, etc.) -the truth value of a property being unknown for those instances too complex to be analyzed. Finally, the form gives information about the size of the reachable marking graph of each instance (number of markings, number of transition firings, maximal number of tokens per place and per marking) -possibly with unknown or approximate answers. The contents of the model forms evolve regularly, as new information and properties are added every year.

Basically, each form is a L A T E X file to be manually written by the person who submitted the model. This approach raises practical issues, as the numerical and Boolean answers provided for sizes and properties are sometimes incomplete or erroneous. Moreover, the annual addition of new models in the collection and of new properties in the forms makes manual handling cumbersome and error-prone. To address this issue, an automated tool chain has been progressively developed by the MCC model team: each P/T net instance is processed in sequence by two tools, PNML2NUPN [START_REF]PNML2NUPN[END_REF] and CAESAR.BDD [1], to produce a L A T E X file containing numerical and Boolean an-swers that help checking and completing manually-written model forms. This tool chain works satisfactorily, although with three limitations: (i) it does not handle colored-net instances -but their structural and behavioural properties, as well as the size of their marking graphs, are expected to be the same as for the corresponding unfolded P/T-net instances; (ii) for P/T nets that are not one-safe, answers given for behavioural properties may be approximate; (iii) for large and/or complex instances, behavioural properties also get approximate answers, as the symbolic (BDD-based) marking-graph exploration (done using the CUDD library [2]) is halted by a timeout after a few minutes.

Nested-Unit Petri Nets. A large proportion (43%) of the MCC models are one-safe P/T nets (also called elementary nets or condition-event systems). Even if such nets are simpler than colored nets or general P/T nets, most related verification problems are computationally hard, namely PSPACE-complete [START_REF] Cheng | Complexity Results for 1-safe Nets[END_REF]; indeed, there are many one-safe P/T-net instances that no participating tool has been able to entirely analyze so far.

A plausible explanation for such difficulties is that verification is made harder because relevant information about the real systems to be modelled has been lost when formalizing these systems as low-level Petri nets: preserving such information all along the modelling and analysis steps would hopefully make verification easier. This is the motivation behind Nested-Unit Petri Nets (NUPN for short) [START_REF] Garavel | Nested-Unit Petri Nets: A Structural Means to Increase Efficiency and Scalability of Verification on Elementary Nets[END_REF], a P/T-net extension that retains structural and hierarchical information by recursively expressing a net in terms of parallel and sequential compositions.

In 2015, the NUPN paradigm was adopted by the MCC team to foster progress in the analysis of large models. A PNML tool-specific extension 3 has been defined to incorporate NUPN information in PNML files, and 8 models (totalling 21 instances) have been enriched with NUPN information. At present, 4 model checkers (CAESAR.BDD, ITS-Tools, LTSmin, and pnmc) are already able to exploit the NUPN information.

Formulas

Tools competing in the MCC are evaluated over four categories of verification tasks: state-space generation, reachability analysis, LTL analysis, and CTL analysis. To maximize tool participation, we further divided the three latter categories into subcategories containing only formulas with a restricted syntax. In 2015, we progressed in the consolidation of these formula languages and provided simplified BNF grammars for each (sub)category, still preserving backward compatibility with previous MCC editions.

Each tool developer may choose in which categories/sub-categories the tool participates. For each model instance and each subcategory, 16 formulas are automatically generated and stored into a single XML file (of which a textual version is also provided for the convenience of tool developers). Each tool participating in the corresponding subcategory is requested to evaluate, on the corresponding instance, all or part of the formulas contained in the XML file.

One of the main criticisms about the formulas generated for MCC'2014 was that many of these formulas were trivially true or false, or could be decided immediately by examining only the initial state or a small fragment of the state space. Generating formulas that are "harder" to model check (not to say realistic verification goals) is a difficult problem, at least because: (i) we need a large amount of formulas (in 2014, we used about 40 000 formulas), and (ii) we usually have little knowledge about the expected properties of the MCC models gathered over the years (namely, only a few models come with associated "meaningful" properties).

To improve the quality of formulas, in 2015 we adopted a new strategy for their generation. Using the grammar of each category, we generated random formulas of up to certain depth (7 operators). We subsequently filtered all generated formulas, in two steps. First, we used SAT solving to filter out formulas being equivalent to true or false independently of the model. Each formula passing the SAT filter was then submitted to SMC [START_REF] Smc | [END_REF], a CTL bounded model checker that we developed ad-hoc for the competition. If SMC was able to decide the satisfiability of the formula by examining only the first 1000 reachable states (computed using BFS exploration), then we discarded the formula. Otherwise, we considered the formula to be hard enough and included it in the XML file. The process continued until either we found 16 hard formulas or we had examined 320 (= 20 × 16) random formulas. In the latter case, we completed the file with (up to 16) random formulas, which might be fairly easy to solve; such a situation mainly happened for small nets having less than 1000 reachable states, or when the formula syntax in a given category made it difficult to find hard formulas -the particular grammar chosen for formula has an impact in this respect, and this is indeed something to work on for the next editions of the MCC.

Among the formulas submitted to SMC, 33.0% were declared satisfiable, 57.4% unsatisfiable, and 9.6% hard. In other words, less than 10% of the generated random formulas were actually retained. About 11.3% of the XML files had to be completed with at least one random formula, leading to a proportion of 8.7% non-hard formulas.

The generation-and-filtering strategy took about two days of CPU to produce the 65 490 formulas used in the 2015 edition. The higher quality of this generation process is confirmed by the increase of CPU time required to process all examinations on all models for all tools. Compared to MCC'2014 (157 days of CPU time), we observed a ten-time increase for MCC'2015 (1541 days of CPU time) where only a two-time increase was expected, based on the new models and larger set of formulas introduced in 2015. Such a larger increase, which actually matches the aforementioned filtering rate (9.6%), clearly indicates that the new strategy produces more demanding formulas that better exercise the capabilities of model checkers.

Monitoring Environment and Experimental Conditions

Due to, at least, the growth of the MCC collection of benchmarks, the number of required executions to evaluate tools is increasing every year. From 54 293 in 2013, it grew to 83 308 in 2014, and to 169 078 in 2015. Such a number of executions thus requires a dedicated software environment that can take benefits from recent multi-core machines and powerful clusters. Moreover, we need to measure key aspects of computation, such as CPU or peak memory consumption, in the least intrusive way.

Since the second edition of the Model Checking Contest in 2012, we have been confident that relying on virtual machines to operate tools is a suitable solution. To achieve this in an automated way, we developed BenchKit [START_REF] Kordon | BenchKit, a Tool for Massive Concurrent Benchmarking[END_REF], a software technology (based on QEMU) for measuring time and memory during tool execution. First used during MCC'2013, BenchKit was then enhanced for MCC'2014 with new management facilities to ease its operation and, for MCC'2015, with the possibility to assign several cores to virtual machines.

BenchKit is operated using configuration files that define consistent sets of runs. Each run represents an execution, i.e., one tool performing one examination on one (instance of a) model. To cope with the ever-increasing need for CPU, we used several machines in 2015, namely: -bluewhale03 (Univ. Geneva), a 40-core, 2.8-GHz machine with 512 GB of RAM, ebro (Univ. Rostock), a 64-core, 2.7-GHz machine with 1024 GB of RAM, quadhexa-2 (Univ. Paris Ouest Nanterre), a 24-core, 2.66-GHz machine with 128 GB of RAM, small (Sorbonne Univ., UPMC), a portion of a cluster consisting of five 24-core, 2.4-GHz nodes with 64 GB of RAM each.

These powerful computing machinery (152 cores in total) enabled each tool to be run in a 64-bit virtual machine with 16 GB of memory and either one core (for sequential tools) or four cores (for parallel tools).

To enable relevant comparison of executions between tools, we divided the set of runs into several consistent subsets. Within each subset, all the examinations related to a given model have been executed on the same machine for all the tools. Then, from the raw data produced by BenchKit (execution time, total CPU, memory peak, sampling of memory and CPU usage), we automatically generated summary charts in HTML format. Fig. 1 provides an example of such a chart comparing two tools in a given examination. Such charts are produced for all the possible tool comparisons and tell how well the verification techniques implemented in the tools scale with the benchmarks under study.

The post-analysis scripts that aggregate data, generate summary HTML pages, and compute scores for the contest are implemented using 15 kLOC of Ada and a bit of bash. BenchKit itself consists in approximatively 1 kLOC of bash.

Participating Tools

Ten tools were submitted in 2015: Cunf [3], GreatSPN-Meddly [4], ITS-Tools [5], LoLA 2.0 [6], LTSmin [7], Marcie [8], pnmc [9], PNXDD [11], StrataGEM 0.5.0 [13], and TAPAAL [14] (four variants, two of them using parallel computing). Table 2 lists the participating tools, and indicates which classes of Petri nets and which verification techniques are supported by each tool: parallel computing, CEGAR, use of decision diagrams (symbolic approaches), explicit model checking, unfolding, use of SAT/SMT solvers, use of state-compression techniques (other than decision diagrams), exploitation of state equations, static variable reordering, use of structural reductions, use of stubborn sets, exploitation of symmetries, use of topological or structural information, unfolding transformation into equivalent P/T nets (for colored nets), and use of Nested-Unit Petri Net (NUPN) information when available.

Results

All examinations on the models were processed in the following way. First, tools were invoked on the "known" models (i.e. those of past years). Then, they were confronted with the "scrambled" models that are the "known" ones but presented as new ones. Finally, they had to process the new "surprise" models (those submitted by the community in 2015). This led to 169 078 runs, from which we extracted a much larger number of computed values (e.g., size of the state space, maximum number of tokens per marking, result of a given formula) that were stored and compiled in order to : (i) evaluate the reliability of the tools, and then (ii) compute an overall score.

Computation of the Reliability of the Tools. Ensuring the reliability of a verification tool is important. We thus decided to use the Model Checking Contest for this purpose too. However, determining if the answers issued by tools are correct is a difficult question, due to the lack of a reference verified result. So, we compare the results of each tool to those of other tools, relying upon a majority-based approach to compute the scores. This is done in three steps: (i) identifying the set of values for which a majority of tools agree; we need at least three agreeing tools to select such a value, (ii) identifying, for each tool, the c computed values that are produced within this set, and (iii) checking, within these c values, the n ones that are correct by the identified majority.

Table 3 (central part) summarizes the reliability rate we computed for each tool. The four columns provide the following data: the computed rate, the number of correct values (C), the number of values in the selected set for the tool (S), and, as an indication, the number of examinations the tool participates in (i.e. provides at least an answer for one instance of a model). Then, the reliability is trivially the ratio C/S. Some outcomes need an explanation. There are tools with a lower reliability rate such as GreatSPN-Meddly, ITS-Tools, LTSmin, or one of the TAPAAL variants. In fact, these tools exhibited bugs introduced in recent development (TAPAAL), or translation mistakes (in formulas for GreatSPN-Meddly and ITS-Tools, in the PNML import for LTSmin). Such bugs could not be identified by their developers before the full benchmark was passed and compared to other tools outputs. We believe this is a valuable information for tool developers, helping them to increase the confidence when implementing new algorithms.

In some cases, we also have tools for which the reliability rate was computed on a very small dataset. Usually, these tools only participate in one examination. For PNXDD, we discovered too late a configuration problem in the virtual machine that was only standing for scrambled and surprise models. We regret this problem was not detected during the qualification procedure. StrataGEM, and pnmc also participated in 3. Summary of results -grey cells mean that the tool was among the top-three winners. When there are several variants for a tool, only the best variant is considered for the podium. only one examination, which probably makes the tool simpler, and thus easier to be more reliable.

We think the procedure is fair, even if rare cases where the majority of tools might be wrong cannot be excluded. Some tool developers have reported that they included the outputs of the MCC in their own built system to increase the reliability of their tool, which will most probably lead, next year, to a better reliability rate.

Computation of the Overall Score for Each Tool. All computed values were then used to elaborate scores. Only correct values were considered and one wrong value disqualifies the tool for the considered examination on the instance of the model. When results were uncertain (i.e., only two distinct answers), the respective reliabilities of the tools were compared to decide which tool was right. We also discarded lonely answers from tools whose reliability was less than 0.9. Then, outputs were weighted: ×1 for "known" models, ×2 for "scrambled" ones and ×3 for "surprise" ones. Finally, points were summed up for each tool and for each examination category: state-space generation, reachability formulas, CTL formulas, and LTL formulas.

The right part of Table 3 summarizes this result. First, no tool did compete on the LTL formulas, which is unfortunate. Marcie won the state-space generation, followed by pnmc and ITS-Tools, which scaled well but exhibited low reliability in many cases. LoLA won the reachability examination, followed by TAPAAL (the "classical" sequential version), and Marcie. For CTL formulas, Marcie is the only participant (and winner) because a bug was signaled too late by the GreatSPN developers in the CTL translation; they preferred to withdraw the tool from this competition.

We are aware that a single table cannot express all details 4 found in the hundreds of thousands charts and tables generated automatically from execution runs. However, we believe such results are useful for the community to: (i) recognize the great effort of building prototype tools and use them for the verification of real systems, (ii) identify the forces and weaknesses of these tools, and (iii) participate in their improvement. In particular, we note that most of the winning tools share a particularity: they combine several techniques, as shown in Table 2.

Conclusion -Future Evolutions

All in one, the fifth edition of the Model Checking Contest was successful and ran smoothly, despite the many novelties introduced in 2015: support of Nested-Unit Petri Nets (to open the contest to verification tools operating on concurrent state machines), new XML syntax for temporal-logic formulas, enhanced techniques to automatically generate non-trivial formulas, and virtual machines with multiple cores in BenchKit. The main disappointment for the MCC team was the absence of submitted tools participating in the LTL category: this could perhaps be interpreted as a sign of decreasing academic interest in linear-time formulas and, if confirmed over the coming years, a pragmatic answer to the longstanding linear-time vs branching-time dilemma.

The following developments are planned for the next editions of the Model Checking Contest:

-Concerning the models: future calls will enrich the collection with new models; the forms will be made more detailed by adding a few missing model properties (e.g., extended free choice); an online repository is being developed, which will offer easy access (using web browsing or programming interfaces) to the collection of models; the proportion of NUPN models is expected to grow, as work is going on to automatically synthesize a NUPN structure for "flat" 1-safe nets. -Concerning the formulas: the community will be actively encouraged to submit properties together with the surprise models, so as to have more "semantically meaningful" formulas associated to each model; the grammar used to describing the syntax of certain categories of formulas will be improved; the generation algorithms will be improved to produce "higher-quality" formulas, even at the cost of discarding some contest categories if necessary. -Concerning the contest itself: new procedures will be added for an early detection of issues similar to that discovered too late with the PNXDD tool (see Sec. 6); also, if tool developers agree, negative scores could be introduced to penalize tools giving incorrect verdicts.

Fig. 1 .

 1 Fig. 1. Example of comparison chart between Tool 1 and Tool 2 . Each point represents an execution and, if below (resp. above) the diagonal, indicates that Tool 1 (resp. Tool 2 ) has won. Faulty executions are displayed on the top and right lines. The grey color means the missing tool did not compete, and the red color means the tool was disagreeing with the majority of other tools.

Table 2 .

 2 List of participating tools in 2015 and the techniques they implement.

	Tool information						Verification techniques
	Name	Institution	P/T nets	Colored nets	Parallel comp.	CEGAR	Dec. Diag.	Explicit	Net Unfold	SAT/SMT	Stat. Compr.	State Eq.	Stat. Var. Re.	Struct. Red.	Stub. Sets.	Symmetries	Topological	Unfold. to P/T	Use NUPN
	Cunf	Univ. Paris 13 (FR)
	GreatSPN-Meddly Univ. Torino (IT)
	ITS-Tools	UPMC (FR)
	LoLA 2.0	Univ. Rostock (DE)
	LTSmin	Univ. Twente (NL)
	Marcie	Univ. Cottbus (DE)
	pnmc	IRT St-Exupery (FR)
	PNXDD	UPMC (FR)
	StrataGEM	Univ. Geneva (CH)
	TAPAAL	Univ. Aalborg (DK)

The collection of benchmarks is available from http://mcc.lip6.fr/models.php

The list of publications is available from http://mcc.lip6.fr

See http://mcc.lip6.fr/nupn.php

See http://mcc.lip6.fr/2015/results.php

Acknowledgments. The MCC team heartily thanks all scientists who proposed new models in 2015 and the tool developers for their participation effort. Acknowledgments are also due to the Universities of Rostock, Geneva, Paris Ouest Nanterre La Défense, and Pierre & Marie Curie for giving access to the powerful computers needed to operate the contest. We are grateful to Hugues Evrard for proof-reading the present paper.


V V Protasov

Nesvizhevsky

G Pignol
email: pignol@lpsc.in2p3.fr

K Protasov

A note on spontaneous emission of gravitons by a quantum bouncer

Spontaneous emission of graviton rates for the quantum bouncer states are evaluated.

The quantum problem of a ball bouncing above a ideal mirror was considered a long time ago, however just as a mere exercise of elementary quantum mechanics. Things did change since the quantization of energy of Ultra Cold Neutrons (UCN) bouncing above a mirror in the Earth's gravitational field had been demonstrated in an experiment performed at the Institute Laue Langevin (ILL) [1,2,3]. However, this effect does not demonstrate any quantum behavior of the gravitational field itself. In analogy with electrodynamics, the observation of spectral lines in atoms shows the quantum behavior of electrons, but does not provide any clue concerning the possible quantization of the electromagnetic field. What does provide a clue is the observation of spontaneous decay of an excited state, for instance, which can only be explained in terms of photon emission.

So the observation of spontaneous decay of an excited state in the ILL experiment would be of interest, since it would be a Planck-scale physics effect. Nevertheless, the decay rate is expected to be low, and the purpose of this letter is to estimate it.

First we will set notations of the quantum bouncer problem, focusing on its physical implementation, that is, neutrons falling on the Earth's gravitational field above a perfect mirror. We then derive the spontaneous emission rate in a semi-classical approach.

The stationary Schrödinger equation for the vertical motion (z axis) of the mass m quantum bouncer is:

- 2 2m d 2 ψ dz 2 + mgz ψ = E ψ. (1)
The boundary condition due to the presence of the perfect mirror at z = 0 is ψ(z = 0) = 0. The characteristic length and the characteristic energy of this problem are:

z 0 = 2 2m 2 g 1/3 = 5.87 µm (2)
E 0 = mgz 0 = 0.60 peV. (3)
The eigenproblem (1) can be solved in terms of the first Airy function Ai(X), which has an infinite number of negative zeros, denoted by {-λ 1 , -λ 2 , . . .} in decreasing order. The energy of the stationary states is equal to E n = E 0 λ n , and the wave function of the n th state is:

ψ n (z) = C n Ai z z 0 -λ n θ(z) (4)
where C n normalizes the probability to find a neutron anywhere to 1. The sequence of the Airy function zeros has no simple analytic expression, but using the Bohr-Sommerfield rules, we find a fairly good approximation of this sequence:

λ n ≈ 3π 8 (4n -1) 2/3 , ( 5
)
this approximation is known to be very good even for the lowest states.

In order to evaluate the rate for a bouncer to make a transition k → n, we will follow the semi-classical analysis as in ref. [START_REF] Weinberg | Gravitation and cosmology[END_REF]. This procedure consists in deriving the classical gravitational power P emitted by an oscillating quadrupole Q cos(ωt), and by replacing the classical quadrupole by the quantum quadrupole moment for the transition k → n which reads

Q kn = m k|ẑ 2 |n , (6)
The quantum mechanical transition rate is:

Γ sp k→n = P ω kn = 4 15
ω 5 kn M 2 P l c 4 Q 2 kn , (7)
M P l is the Planck mass and ω kn = (E k -E n )/ is the angular frequency of the transition. This formula is valid if the quadrupole approximation ω kn z k ≪ c holds, i.e. k ≪ 10 8 . This semi-classical derivation for spontaneous emission of gravitons was recently considered in ref. [START_REF] Boughn | [END_REF] for atomic hydrogen. The authors provide a field theory derivation of Eq. ( 7) and show that it satisfies the detailed balance.

Since the quadrupole matrix elements for the quantum bouncer are known explicitly in terms of the Airy function zeros [6], we get:

k|ẑ 2 |n = 24(-1) k-n+1 (λ k -λ n ) 4 z 2 0 . (8)
Or, the transition probability is equal to

Γ sp k→n = 512 5 1 (λ k -λ n ) 3 m M P l 2 E 5 0 z 4 0 c ( c) 5 = 5 × 10 -77 s -1 (λ k -λ n ) 3 . ( 9
)
For the two lowest quantum states λ 2 -λ 1 = 1, 75 and the probability of the spontaneous graviton emission is as low as:

Γ sp 2→1 ∼ 10 -77 s -1 . ( 10
)
This effect is unfortunately experimentally too rare to observe.