The IAS-MEEG Package: A Flexible Inverse Source Reconstruction Platform for Reconstruction and Visualization of Brain Activity from M/EEG Data

We present a standalone Matlab software platform complete with visualization for the reconstruction of the neural activity in the brain from MEG or EEG data. The underlying inversion combines hierarchical Bayesian models and Krylov subspace iterative least squares solvers. The Bayesian framework of the underlying inversion algorithm allows to account for anatomical information and possible a priori belief about the focality of the reconstruction. The computational efficiency makes the software suitable for the reconstruction of lengthy time series on standard computing equipment. The algorithm requires minimal user provided input parameters, although the user can express the desired focality and accuracy of the solution. The code has been designed so as to favor the parallelization performed automatically by Matlab, according to the resources of the host computer. We demonstrate the flexibility of the platform by reconstructing activity patterns with supports of different sizes from MEG and EEG data. Moreover, we show that the software reconstructs well activity patches located either in the subcortical brain structures or on the cortex. The inverse solver and visualization modules can be used either individually or in combination. We also provide a version of the inverse solver that can be used within Brainstorm toolbox. All the software is available online by Github, including the Brainstorm plugin, with accompanying documentation and test data.


Introduction
The recording of the extracranial magnetic field, or alternatively, the voltage distribution on the scalp, induced by the concerted firing of bundles of neuron in the brain constitute the data of MagnetoEnchephaloGraphy (MEG) or Electro-EncephaloGraphy (EEG), two modalities for monitoring brain activity in a totally non-invasive manner with an exquisite time resolution in the millisecond range. It is because of its high temporal resolution that the M/EEG modalities are the tools of choice to investigate and localize brain phenomena that occur within a short time interval, such as the initiation of an epileptic seizure.
The challenges in the use of the M/EEG modalities come from the weakness of the signal, and ill-posedness of the inverse problem of mapping the boundary data to the brain region. In particular, the number of M/EEG channels is of the order of one to a few hundreds, while the brain activity is typically represented in the form of an ensemble of current dipoles at a number of locations in the brain ranging in tens of thousands, distributed on the cortex as well as in the internal structures. The difficulties associated with the solution of the M/EEG inverse problem include the low signal-to-noise ratio, the high sensitivity of the solution to the perturbations in the data, and the non-uniqueness of the solution (Hämäläinen et al. 1993;Baillet et al. 2001;Brette and Destexhe 2012), thus requiring that the data is augmented with additional information in the form of a regularization term or a prior. An additional byproduct of the non-uniqueness of the solution is that most classical inversion methods tend to favor sources closer to the sensors, thus giving preference to cortical activity over activation in the deeper brain.
Inverse problems can be recast in the form of Bayesian inference problems, and several algorithms based on Bayesian hierarchical models can be found in the literature (Auranen et al. 2005;Calvetti et al. 2009;Henson et al. 2009Henson et al. , 2010Kiebel et al. 2008;Lucka et al. 2012;Lopez et al. 2014;Mattout et al. 2006;Nummenmaa et al. 2007a, b;Owen et al. 2012;Sato et al. 2004;Stephan et al. 2009;Trujillo-Barreto et al. 2004;Wipf and Nagarajan 2009;Wipf et al. 2010). The inverse solver algorithm discussed in this paper (Calvetti et al. 2009(Calvetti et al. , 2015 is based on a model where the unknown dipoles are conditionally independent Gaussian random variables whose variances, also unknown, are modeled as random variables following a gamma distribution. The hierarchical structure of the prior, particularly suitable for modeling focal activity patterns, effectively doubles the number of unknowns, as now both the dipoles and their variances need to be estimated. The negative logarithm of the resulting posterior density, or Gibbs energy, has been shown to be globally convex (Calvetti et al. 2009), hence to have a unique minimizer. The Iterative Alternating Sequential (IAS) algorithm computes an approximation of the Maximum a Posteriori (MAP) single estimate of the posterior. Each iteration of the algorithm consists of a sequence of two steps, one requiring the solution of a least squares problem to update the estimate of the dipoles and the other consisting of a formula evaluation to update the estimate of their variances.
In addition to the theoretical advantage of guaranteed unique MAP estimate, the computational advantage of the software lies in the implementation of the optimization algorithm. The two alternating steps consist of a closed formula update of the variances, and a linear least squares problem for updating the dipoles. The latter step, that is potentially time consuming, is implemented in IAS by using a Krylov subspace iterative method based on the Conjugate Gradient for Least Squares (CGLS) with a particular right preconditioner coming from a symmetric factorization of the precision matrix of the prior, therefore referred to as priorconditioner, and a left preconditioner coming from a symmetric factorization of the precision matrix of the noise. Unlike in the traditional use of iterative solvers for linear systems where the prior provides a regularization, the iterative algorithm in IAS solves a reduced low rank linear system informed by the prior, and it is equipped with a suitable early stopping rule, acting as a regularization. This novel approach, the priorconditioned Krylov solver scheme, reduces dramatically the computational cost and has the advantage of converging very fast, typically requiring only few (in the order of tens) matrix-vector products with the lead field matrix, is more than an alternative way of solving a linear system, since the solution with early stopping depends nonlinearly on the data (Calvetti et al. 2018). The subject-specific anatomical information based on the segmented MRI images is encoded in the right preconditioner. Anatomical information about the brain to augment the data has been previously used in MEG algorithms, see, e.g., in .
The popularity of standard public domain M/EEG solvers, e.g., Minimum Norm Estimate (MNE) (Hämäläinen and Ilmoniemi 1984;) and LORETA (Pascual-Marqui 1999) is partly due to the fact that they can be used with minimal user intervention. This observation has motivated our effort to reduce the number of user-provided parameters in the present IAS-MEEG to a minimum, with a clear physical interpretation. The key parameter is an estimate of the signal-to-noise ratio (SNR). Another optional input parameter controls the focality of the solution. In Calvetti et al. (2009), the sparsity-promoting propensity of the IAS algorithm, and its relation to some other algorithms favoring sparse solutions (Gorodnitsky and Rao 1997;Nagarajan et al. 2006;Uutela et al. 1999) was highlighted. In particular, it was shown that the shape paraof the gamma hyperprior distribution controls the sparsity, and, furthermore, in , it was shown that at an appropriate limit of the shape parameter, the IAS solution converges towards the weighted Minimum Current Estimate (MCE) (Uutela et al. 1999), providing a good intuition for the role of the shape parameter.
Sensitivity weighting, or depth weighting, has been commonly used in connection with, e.g., MNE and MCE methods to overcome the bias of minimization based methods towards sources closer to the measuring devices. Recently, a proper Bayesian interpretation of the sensitivity weighting has been provided in . In the present version of the algorithm, the sensitivity weighting is automatic, arising from a very natural exchangeability argument , requiring that for the given SNR, all dipole configurations with equal cardinality of the support should a priori be equally probable.
In , the performance of the IAS-MEEG package for the identification of an active brain region ranging from the cortex to the deep structures of the basal ganglia was systematically compared to that of three other standard inversion methods provided in Brainstorm (Tadel et al. 2011): wMNE , dSPM (Dale et al. 2000), and sLORETA (Pascual-Marqui 1999;Wagner et al. 2004). Extensive statistical evaluations demonstrated that IAS-MEEG is very well suited for the localization of focal brain activity in both cortical and subcortical regions, and its ability to identify an activation pattern even in the presence of disturbances due to internal brain noise was confirmed.
In this paper we present a standalone software platform for the reconstruction of the neural activity in the brain from MEG or EEG data based on the IAS algorithm, and a visualization of the result as sagittal, coronal and axial slices. The 1 3 reduction in computing time due to the iterative priorconditioned Krylov subspace least squares solver makes the software particularly suitable for the reconstruction of lengthy time series on standard computing equipment. The algorithm requires minimal user provided input parameters, effectively reduced to an estimate of the signal-to-noise ratio. Optionally, a parameter controlling the focality of the solution, the desired accuracy of the solution, and a parameter related to expected maximal strength of the sources may be entered. The input data consist of the MRI and M/EEG measurements, cleaned of obvious artifacts due, e.g., to heartbeat or eye motion, a discretized source space, and the lead field matrix. The computational procedure consists of the following processing stages: (1) Set up the anatomical prior from information in the MRI data; (2) Initialize the vector of the variances of the hyperpriors for the dipoles according the Bayesian sensitivity weighting formula; (3) Initialize the variances of the dipoles; (4) Update alternatively the estimate of the current dipole moments and their variances until convergence to a specified accuracy has been achieved; (5) Visualize the reconstructed activity.
One of the novelties of the current article is the integration of results from several earlier works, making the model scaling and parameter selection automatic, rendering the necessary user interference minimal, limited to the input of an estimate of the signal-to-noise ratio. However, the algorithm leaves the option to control some of the features such as the parameter controlling the required source focality in the solution. Another novel aspect is the integration of EEG and MEG modalities in the same IAS algorithm, which makes the platform particularly flexible.
The paper is organized as follows. After a brief review of the M/EEG inverse problem from a Bayesian perspective in Section "Materials and Methods", we recall the main step of IAS algorithm in "The IAS Algorithm" section and discuss the physiological meaning of its parameters. Section 4 describes the lightweight visualization tool that we provide together with the algorithm, and in Section 5 the scripts contained in the IAS-MEEG package are described. Finally, in Section 6 we present some computed examples with synthetic and real data.

Symbols and Notation
For the convenience of the reader, we begin by stating the assumptions and then establish the notation that we will use in the following sections.
The source space is assumed to have n dipoles, some of which have a natural preferred orientation determined by the underlying neuroanatomy. For the j-th dipole, 1 ≤ j ≤ n , we denote: When needed, the dipole variables are collected in the vectors: The local anatomical prior matrices are used to construct the block diagonal matrix that admits a Cholesky factorization as where is an upper triangular matrix. In the following we will also use the matrix and its factorization In the solution of the M/EEG inverse problem we use the notation We will denote a time slice of data by the vector b ∈ ℝ m , with m the number of M/EEG sensors.
Finally, for a vector V ∈ ℝ n and a matrix ∈ ℝ n×n we use the norms: and r j ∈ ℝ 3 point coordinates, ∈ ℝ m×3n the lead field matrix, B ∈ ℝ m×t a dataset of t time slices, ∈ ℝ m the noise vector, ∈ ℝ m×m the noise covariance matrix.

The Bayesian Hierarchical Model
The cortical and sub-cortical surfaces of the brain are extracted from the MR images of the subject under study and discretized in a regular triangulation whose nodes form the source space. A current dipole is located in each point r j , 1 ≤ j ≤ n , of the source space. Let q j , 1 ≤ j ≤ n , be the moment of the j-th dipole. The primary unknown is the current dipole moment vector Q = [q 1 , … , q n ] . Each dipole has a preferred direction j that can be extracted from the MRI. The local prior variance of the amplitude of each dipole j is modeled as a random variable following the Gamma distribution with hyperparameters ( j , * j ): The observation vector b ∈ ℝ m and the dipole moment vector Q ∈ ℝ 3n are assumed to be linearly related: where ∈ ℝ m×3n is the lead field matrix and ∈ ℝ m is the observation noise vector. Modeling the noise term as a zero mean Gaussian random variable, ∼ N(0, ) , where ∈ ℝ m×m is the noise covariance matrix, the likelihood density of b conditional on Q can be written as We define the Bayesian hierarchical prior model of the activity in position r j as where j ∈ ℝ 3×3 is the local anatomical prior matrix and ( j , j , j ) is a local orthonormal frame at r j .
Assuming a priori that the dipoles are conditionally independent leads to a conditionally Gaussian prior model with the hyperprior model It follows from Bayes' theorem that the posterior distribution is of the form The Maximum A Posteriori (MAP) estimate of both the dipole moment vector Q and corresponding variance vector is obtained minimizing the energy function by the IAS algorithm. In fact, the minimization with respect to Q affects only part (a) and reduces to a quadratic minimization problem that can be solved very efficiently using a priorconditioned CGLS algorithm, while the minimization with respect to the hyperparameters j affects only part (b) and admits a solution in closed form.

The IAS Algorithm
The Algorithm The IAS algorithm minimizes the energy function E(Q, ) by proceeding through a sequence of iteration steps alternating the minimization with respect to Q and the minimization with respect to . The algorithm is initialized by assigning to a given value, for instance, * . The first step is to minimize part (a) of the energy function (2) with fixed = * . The MAP estimate of Q given = * is obtained by solving the minimization problem The minimizer can be found by solving an associated linear system in the least squares sense, and an efficient implementation can be done using a priorconditioned CGLS algorithm, see (Calvetti et al. 2015) for details. Once Q est has been computed, the updated value of is obtained minimizing part (b) of (2) keeping Q = Q est fixed, (3) This problem has the analytical solution where j = j − 5∕2 . Subsequently, a new estimate of Q is obtained by solving problem (3) again with = est . Thus, the minimization problems (3) and (4) are solved alternatingly until a convergence criterion is met. The flowchart of the IAS algorithm is shown in the left of Fig. 1 with the pseudocode of the algorithm shown in the right.

Choice of the Parameters
While the IAS algorithm depends on two parameter vectors, the shape parameter vector and the scaling parameter vector * , the user needs only to enter an estimate for the signal-to-noise ratio. Internally to the algorithm, the scaling parameter is automatically computed, and the shape parameter is in practice reduced to a single scalar input that can be set at a default value. The description is given below. The shape parameters j > 5∕2 , 1 ≤ j ≤ n , control the focality of the reconstructed sources: values close to the lower bound favor sparse solutions while larger values favor distributed solutions. For the sake of simplicity, we assume that the same value of j = is used for all dipoles in the source space. Recalling that j = j − 5∕2 , 1 ≤ j ≤ n , we set j = , 1 ≤ j ≤ n . Reasonable values of are in the interval [0.0001, 0.1]: when = 0.0001 the algorithm favors current density estimates with narrow support, while for = 0.1 the support of the estimated density is wider.
The scaling parameters * j , 1 ≤ j ≤ n , are related to the expected value of the variance j of the j-th dipole. By using a Bayesian argument that is based on requiring exchangeability of dipoles to explain the signal-to-noise ratio shows that * j needs to be chosen inversely proportional to the sensitivity with respect to the j-th dipole. Its value can be evaluated explicitly by the formula where P is the power of the signal, is the noise covariance matrix, j ∈ ℝ m×3 is the local lead field matrix and ‖ ⋅ ‖ F est = argmin denotes the Frobenius norm . Recalling the definition of the signal-to-noise ratio (SNR) as we arrive at the formula Assign the parameters η = β − 5/2, Θ * , and τ Fig. 1 The IAS flowchart (left) and the IAS algorithm (right). The outer loop on k is performed through statements 6-25; the inner loop (statements 12-20) solves the minimization problem for Q k+1 by the preconditioned CGLS algorithm; the estimated Q k+1 is computed by statement 21; the estimated k+1 is computed by statements 22-24. Observe that while the IAS algorithm requires the vector * as input, the user specifies only the signal-to-noise ratio, and * is automatically generated based on the SNR-exchangeability argument. The parameters and have default values, and specifying them by the user is optional In , the formula was related to the estimate of active focal sources via the exchangeability argument.

Time Series
In case of time dependent data, the IAS algorithm is applied to each individual time slice. The time slices can be processed individually, or treated as a time series. In the latter case, the hyperparameter vector is initialized to * in the first time step, and set equal to the final estimate in the previous time step in each subsequent time steps. The rationale for the dynamic initialization of is that since we do not expect significant changes in the current density from one time instance to the next, it is reasonable to assume the variance not to vary much either, thus making the variance at time step t a good estimate of the variance at time step t + 1 . Computed examples show that in this way, the outer iteration loop typically requires only a few iterations to reach convergence.

Visualization of the Activity Map
To provide a standalone platform, the software is equipped with a computationally lightweight visualization tool based on basic graphic packages of Matlab. The visualization provides a sliced view of the reconstructed activity in three standard orthogonal anatomical planes, axial, coronal, and sagittal views. Given the source space in terms of the projections onto the three principal directions, where the default right handed Cartesian coordinates axes are assumed to be in the order x = right , y = front and z = crown (see Fig. 2). The sliced visualization algorithm subdivides the source space in the chosen direction into ten layers. More precisely, for instance in the axial projections, the vertices are parceled as S = [r 1 , r 2 , … , r n ], r j = (x j , y j , z j ), where z = z min + (z max − z min )∕10 with = 0, 1, … , 10 , and z min and z max define the extremal values of the source space components z j . Given the IAS-MEEG reconstruction of the brain activity, the visualization algorithm shows the intensities a j = ‖q j ‖ of the dipoles in each parcel, r j ∈ S axial , by plotting the projection (x j , y j ) as a dot in two dimensions, encoding the intensity in the form of a color map.
In addition to the three sliced views of the full brain, the visualization software allows the user to select a point inside the brain, and produce an activity plot in the form of the single axial, coronal, and sagittal slices containing the selected point. This view is particularly useful if synthetic data corresponding to a focal source are used to monitor how well the algorithm is able to localize a focal source.
The package is coded in Matlab using only basic commands and does not require any of the Matlab toolboxes.
The package comprises four functions: 1. BuildAnatomicalPrior: This function generates the anatomical prior that favors dipoles in the preferred direction. The inputs are: the coordinates of the points in the source space; the normal vectors extracted from the MRI. The user can choose the value of , the relative variance of the components of the dipoles. The default value is 0.05. The output is the matrix ∈ ℝ 3n×3n , that is the Cholesky factor of the anatomical prior covariance matrix. 2. SetParameters: This function scales the lead field matrix and the data, adjusts the truncation of the sensitivities, and returns the scaling vector * together with an estimate for the standard deviation of the noise for whitening. The inputs are: the lead field matrix; the Cholesky factor of the anatomical prior covariance as computed by BuildAnatomicalPrior; a small set of the data; the estimated signal-to-noise ratio; the percentage of the highest * j values to be removed (optional: by default no truncation). Removed values can be checked by graphical inspection asking the function to produce plots. The outputs are: the * vector; the cut-off value used to remove the higher values of * ; the standard deviation of the scaled noise; the scaling factor for the lead field matrix; the scaling factor for the data. _ : This function solves the M/EEG inverse problem by the IAS algorithm described in "The Algorithm section. The inputs are: the lead field matrix; the dataset; the Cholesky factor of the anatomical prior covariance matrix, evaluated in BuildAnatomi-calPrior; the parameters, evaluated in SetParameters; the value of for selecting the focality of the reconstructed activity (default value: 0.01; choose 0.001 for focal sources). The user can choose the maximum number of iterations in the outer and inner loops, n outer and max it , respectively, and the tolerance used in the stopping criterion for . Default values are: n outer =30, max_it=120 and = 0.01. The output is the estimated dipole moment vector for each point of the source space, and a diagnostics matrix indicating the number of inner iterations for each outer iteration, and the relative change in when tolerance is reached.

4.
_ : This function plots the activity map, that is the estimated current intensity of each dipoles (cf. "Visualization of the Activity Map" section). The inputs are: the coordinates of the points in the source space; the estimated intensity vector. The outputs are the plots of the dipole intensity for different sections (axial, coronal and sagittal sections).
In addition, a visualization algorithm is included, which is particularly useful when simulated data are used.

4'
_ : This function plots three sections (axial, coronal and sagittal sections) of the activity map passing through a given point in space. The inputs are: the coordinates of the points in the source space; the estimated intensity vector; the point at which the three views intersect. Optional inputs include the coordinate system specification and the type of marker indicating the intersection of the three views in the plot.
To solve the M/EEG inverse problem by the IAS-MEEG algorithm the four functions must be called sequentially as shown in the flowchart in Fig. 3. The first box of the flowchart shows all the input necessary to run the IAS algorithm. The lead field matrix , the source space and the corresponding normal orientation that need to be given in input can be computed using an available software package (e.g. Brainstorm or Fieldtrip). These inputs are passed as Matlab variables to the different functions as specified in the correspondig box of the flowchart in Fig. 3. We also assume that the input data matrix B contains M/EEG data that were already preprocessed (e.g. filtering, artifact removal) and is also passed as Matlab variable.

Brainstorm Plugin
The IAS algorithm is also available as a Brainstorm plugin, allowing its integration into the Brainstorm workflow. The added value of this integration is that it makes it possible for a user to perform all analysis steps, such as data preprocessing and visualization, within the same toolbox. The three main functions of IAS code, BuildAnatomicalPrior, SetParameters and _ are integrated into the single function _ available in the IAS-MEEG package. To be used within Brainstorm this function must be copied in the Brainstorm folder brainstorm3/toolbox/process/functions and launched as a process. Figure 4 shows the main window for running IAS within Brainstorm with the input parameters to be provided by the user.

MEG Simulated Data
An example of the impact of the choice of shape parameter on the reconstruction is shown in Fig. 5. Here, we used simulated MEG data generated by a patch of activity located in the occipital region of the left hemisphere. The same anatomical information and MEG sensors used for real dataset (cf. "MEG Real data" section) were used. A biological noise, simulated as in , was added to the signal, yielding SNR = 15 . It is clear that when = 0.0001 (see Fig. 5, bottom) the IAS algorithm favors sparse current density estimates respect to a choice of = 0.1 (see Fig. 5, middle) where the estimated density has a more spread distribution.

MEG Real Data
In the documentation (https:// github. com/ IAS-code/ IAS-MEEG/ blob/ master/ IAS_ Demo.m) we provide an example of reconstruction obtained by applying IAS-MEEG algorithm to a real dataset. We used the MEG sample data provided by MNE software package (Gramfort et al. 2014) which includes also the MRI reconstructions created with FreeSurfer. These MEG data were acquired with the Neuromag Vectorview system at MGH/HMS/MIT Athinoula A. Martinos Center Biomedical Imaging. In the protocol experiment, checkerboard patterns were presented into the left and right visual field, interspersed by tones to the left or right ear. The interval between the stimuli was 750ms. Occasionally, a smiley face was presented at the center of the visual field. The subject was asked to press a key with the right index finger as soon as possible after the appearance of the face. In our example we only consider the trials corresponding to the left visual stimulus and perform the averaging on these trials. The MRI data were imported in Brainstorm (Tadel et al. 2011) to generate a source space including both the cortical surface and substructure regions. Following , we adopt for the source space the DBA (Deep Brain Activity) model proposed by Attal et al. (2012), Attal and Schwartz (2013) where the deep regions are modeled either with surfaces or volumes depending on anatomical information. The source space obtained in this manner consists of around 19000 nodes and the lead field matrix was computed using the single layer model implemented in OpenMEEG (Gramfort et al. 2010) software provided by Brainstorm. We reconstructed the neural activity by IAS-MEEG using the following input parameters: SNR = 9 , = 0.01 and cut_off = 0.9 . Figure 6 shows the reconstructed activity for the left visual stimulus at 92ms where evidently an activation in the right occipital region appears. The visual stimulus elicits activity in different regions of the visual cortex, thus to get an overview in time of the different activated regions at different time inputs: R, q 1 (time) , · · · , q n (time) outputs: plots of the activity map Fig. 4 Main window of IAS process within Brainstorm toolbox points, the dipoles with maximum activity were selected and visualized. Figure 7 shows these reconstructed dipoles with the corresponding time trace.

EEG Real Data
To show an example on how to use the IAS-MEEG algorithm within Brainstorm we run the plugin on the EEG sample dataset provided in the Tutorial EEG and Epilepsy https:// neuro image. usc. edu/ brain storm/ Tutor ials/ Epile psy. The tutorial dataset was acquired in a patient who suffered from focal epilepsy at the Epilepsy Center Freiburg, Germany. The EEG data was recorded at 256Hz, using a Neurofile NT digital video-EEG system with 128 channels and a 16-bit A/D converter. The spikes were marked with Brainstorm by the epileptologists at the Epilepsy Center in Freiburg. The individual MRI data were imported in Brainstorm and a cortical source space with around 15000 nodes was created. The lead field matrix was computed using a three layer model (scalp, skull and bran) by OpenMEEG software provided in Brainstorm. The IAS-MEEG algorithm was applied to the EEG data averaged on the marked spikes. Figure 8 shows the IAS localization on epileptic spikes when the following input parameters were used: SNR = 9 , = 0.001 and cut_off = 0.9 . In particular, we used a value of that promotes focality.

Conclusions
In this article we presented the IAS-MEEG package, a standalone, Matlab-based, freely downloadable software for the reconstruction of the neural activity from M/EEG data, with a plugin to integrate it in Brainstorm. The IAS-MEEG package is based on an Iterative Alternating Sequential (IAS) inversion algorithm that combines hierarchical Bayesian models and Krylov subspace iterative least squares solvers. The package is available via GitHub at https:// github. com/ IAS-code/ IAS-MEEG and distributed under a Berkeley Software Distribution (BSD) license. An online documentation is also provided at https:// ias-code. github. io/ IAS-MEEG/ index. html. All core routines are written in standard Matlab language and do not rely on any special packages. Furthermore, the inverse solver and visualization modules can be used either individually or in combination, allowing an easy integration in Brainstorm (Tadel et al. 2011).

Appendix
In this appendix, we review the scaling of the fields and matrices as implemented in the program SetParameters. To scale the lead field matrix, we write the action of the lead Hence, the lead field vector p (j) k indicates how much a unit dipole at position r k contributes to the jth channel. We define a scaling factor indicating the average effect over all dipoles on the channel that is most sensitive to a given dipole. To equilibrate numerically the lead field matrix, we write the forward model (1), indicating the time dependency with a subscript t, 1 ≤ t ≤ T , as To scale the data, we define the mean amplitude over the time series, and scale the Eq. (5) to read or, using the evident notation, The standard deviation sc of the scaled noise is estimated based on the signal-to-noise ratio, and the Eq. (6) is whitened by dividing it by sc .