Computational Intelligence and Neuroscience

Volume 2016 (2016), Article ID 3979547, 11 pages

http://dx.doi.org/10.1155/2016/3979547

## MEG Connectivity and Power Detections with Minimum Norm Estimates Require Different Regularization Parameters

^{1}Psychology Department, University of Montreal, Montreal, QC, Canada H2V 2S9^{2}Department of Computer Science, Pontificia Universidad Católica de Chile, 7820436 Santiago de Chile, Chile^{3}School of Psychology and Interdisciplinary Center for Neurosciences, Pontificia Universidad Católica de Chile, 7820436 Santiago de Chile, Chile^{4}Lyon Neuroscience Research Center, DyCog Team, Inserm U1028, CNRS UMR5292, 69675 Bron Cedex, France^{5}Department of Neuroscience and Biomedical Engineering, Aalto University, 02150 Espoo, Finland^{6}MEG Center, CERMEP, 69675 Bron Cedex, France

Received 15 October 2015; Revised 19 January 2016; Accepted 14 February 2016

Academic Editor: Jussi Tohka

Copyright © 2016 Ana-Sofía Hincapié et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Minimum Norm Estimation (MNE) is an inverse solution method widely used to reconstruct the source time series that underlie magnetoencephalography (MEG) data. MNE addresses the ill-posed nature of MEG source estimation through regularization (e.g., Tikhonov regularization). Selecting the best regularization parameter is a critical step. Generally, once set, it is common practice to keep the same coefficient throughout a study. However, it is yet to be known whether the optimal lambda for spectral power analysis of MEG source data coincides with the optimal regularization for source-level oscillatory coupling analysis. We addressed this question via extensive Monte-Carlo simulations of MEG data, where we generated 21,600 configurations of pairs of coupled sources with varying sizes, signal-to-noise ratio (SNR), and coupling strengths. Then, we searched for the Tikhonov regularization coefficients (lambda) that maximize detection performance for (a) power and (b) coherence. For coherence, the optimal lambda was two orders of magnitude smaller than the best lambda for power. Moreover, we found that the spatial extent of the interacting sources and SNR, but not the extent of coupling, were the main parameters affecting the best choice for lambda. Our findings suggest using less regularization when measuring oscillatory coupling compared to power estimation.

#### 1. Introduction

Healthy brain function is largely mediated by coordinated interactions between neural assemblies in different cortical and subcortical structures. As a result, the study of neuronal processes requires techniques that can reliably measure the spatiotemporal dynamics of large-scale networks. To this end, it is important to assess the modulations of local activations as well as long-range coupling between brain areas. Over recent years, the quantification of neuronal interactions in a given behavioral task or brain state has been the focus of a large body of research, and various techniques are currently used to detect and probe the role of functional connectivity (e.g., [1, 2]). In particular, a widely used technique for the detection of large-scale interactions among neural assemblies is magnetoencephalography (MEG) [3–5]. This is primarily due to its high temporal resolution, which is in the same order of magnitude as the neuronal processes themselves (milliseconds). Unlike electroencephalography (EEG), MEG measures magnetic fields that are less affected by the skull and brain tissue and provides whole-head coverage, which is mandatory for the assessment of large-scale brain networks. The use of MEG has advanced our comprehension of the mechanisms underlying functional brain connectivity [6] involved in sensory, motor, and higher-order cognitive tasks [3, 7–11] higher-order cognitive tasks resting state [6, 12–18].

The ability to measure brain interactions at the source level, rather than between sensor channels, is an important prerequisite if we want to make useful inferences about the anatomical-functional properties of the network. To this end, source reconstruction techniques are applied in order to estimate the spatiotemporal activity of the cortical generators underlying the recorded sensor-level MEG data. Solving this relationship is an ill-posed inverse problem for which numerous methods have been developed [4]. The differences between the various techniques mainly stem from the assumptions they make about the properties of the neural sources and from the way they incorporate various forms of a priori information, if any is available. Conceptually, solving the MEG inverse problem boils down to solving a system of equations that is underdetermined (i.e., no unique solution) since we generally have far more sources (thousands) than measurements (hundreds). Identifying a solution to such a problem can be achieved by imposing constraints on the sources. A typical constraint is to minimize the source power. Among these methods, the Minimum Norm Estimate (MNE) relies on minimizing the L2-norm and is one of the most widely used techniques [4, 7, 8, 18–37]. By contrast, estimates obtained by minimizing the L1-norm are referred to as Minimum-Current Estimates (MCE) [34, 38]. While the L2-norm assumes a Gaussian a priori current distribution, the L1-norm assumes a Laplacian distribution [39].

In principle, MNE looks for a distribution of sources with the minimum (L2-norm) current that can give the best account of the measured data. As the problem is ill-posed, MNE generally uses a regularization procedure that sets the balance between fitting the measured data (minimizing the residual) and minimizing the contributions of noise [22, 40, 41]. The Tikhonov or Wiener regularization [42, 43] and SVD truncation are among the most widely used regularization procedures. Regularization may be considered a necessary evil: it is required to stabilize the solution of the inverse problem, yet too much regularization leads to overly smooth solutions (spatial smearing). Since we do not have a precise model of brain activity and noise, the choice of the optimal amount of regularization is a nontrivial step for which no magical recipe exists [44]. The relationship between optimal regularization and the patterns of underlying generators is still poorly understood. In particular, the effect of regularization on the detection accuracy of Minimum Norm Estimates of source power and interareal source coupling has not been thoroughly investigated. This raises the question of whether one should use the same regularization coefficient for MNE-based power and connectivity analyses, as is generally done, or would one benefit from optimizing the regularization coefficients separately for each analysis? Furthermore, how does the optimal regularization coefficient, in such configurations, depend on sensor-level SNR, source size, or coupling strength?

With these questions in mind, we performed extensive Monte-Carlo simulations, creating over 20,000 pairs of coupled oscillatory time series in MEG source spaces, and computed the resulting surface MEG recordings. Then, we estimated source power and coherence using the MNE framework with variable degrees of Tikhonov regularization. Moreover, we used an approach based on an area under the curve (AUC) to identify the optimal regularization coefficient (lambda) in the case of power and coherence analyses. We found a systematic difference between the optimal lambdas in each analysis: for source-level coherence analysis the optimal lambda was two orders of magnitude smaller than the best lambda for power detection. Lastly, our findings are broken down as a function of SNR, cortical patch size, and corticocortical coupling strength.

#### 2. Methods and Materials

In this section, we first present the MEG inverse problem formulation and the MNE framework. Then, we describe the simulation, reconstruction, and performance assessment procedures.

##### 2.1. The MEG Inverse Problem

The inverse problem looks for an estimation of the active sources () that generates the measurements () recorded at the sensors. The dimension of the columns of accounts for the three components of the source in the , , and directions. According to anatomical observations, the main generators of MEG are located in the grey matter and their orientation is perpendicular to the cortical sheet [45]. Here, we used a constrained orientation approach (perpendicular to the cortical surface), and hence the dimensions of are reduced to . Assuming a linear relationship between the measurements and the active sources, the problem is modeled aswhere is the lead field matrix () and is additive noise applied at the MEG channels (). The lead field matrix describes how each source contributes to the measurements at each sensor, given a specific head conductivity model and a source space. As the number of sources is usually much higher than the number of sensors, the lead field matrix is highly underdetermined and thus not invertible. The estimation of the activity of the sources requires the definition of an inverse operator :where represents the estimated sources () and the superscript denotes matrix transpose.

##### 2.2. Minimum Norm Estimate (MNE)

As the MEG inverse problem is ill-posed, a regularization scheme is needed [22], and one of the most common options is the Tikhonov regularization [42, 43]. MNE calculates an inverse operator by searching for a distribution of sources with minimum currents (L2-norm) that produces an estimation of the measurements () most consistent with the measured data (). The solution is a trade-off between the norm of the estimated regularized sources current and the norm of the quality of the fit they provide to the measurements . Assuming the noise and the sources strength to be normally distributed with zero mean and covariance matrix and , respectively, a general form of the MNE inverse solution can thus be given as [35, 46, 47]where is the Tikhonov regularization parameter. Thus, the inverse operator is defined aswhere the superscript denotes matrix inverse.

The Minimum Norm Estimate embodies the assumption of independently and identically distributed (IID) sources, which corresponds to an identity matrix in the above formula. Alternatively, can incorporate more informed (spatial) assumptions, yielding a so-called weighted Minimum Norm Estimate [28, 35, 44]. Here, we use the general and classical minimum norm solution. The noise covariance matrix was computed from the actual noise which was added to the sensors in each simulation.

##### 2.3. Simulations

We simulated 117s of oscillatory activity in pairs of cortical sources with different degrees of coupling in the alpha band (9–14 Hz). We computed the resulting sensor-level data through forward modeling based on a 275-channel CTF MEG system configuration. Our sources consisted of current dipoles placed at the vertices of a tessellated MNI-Colin 27 cortical surface, which was segmented and tessellated using FreeSurfer [48] and downsampled to 15028 vertices. Different strengths of coupling were obtained by forcing the time series of the second source to have a certain level of coherence with the first source time series. The magnetic fields at the sensors were calculated using a single sphere head model and constraining the orientations of the sources to be normal to the cortical surface. The simulated data were generated using a combination of custom MATLAB code, in addition to functions from Brainstorm [49] and FieldTrip [50] toolboxes. The source time courses were simulated by first setting the base frequency of the oscillator (e.g., 12 Hz) and inducing a small jitter to its instantaneous frequency across time points. The frequency modulated time courses were then generated using an exponential function. This procedure causes fluctuations in the phase relationship between two oscillators with the same base frequency, allowing us to achieve coherence levels below 1. The frequency jittering was performed randomly in a loop until the desired coherence level (e.g., 0.4) between the time courses of the two oscillators was reached.

We randomly selected two locations (seeds) for each simulated pair of sources (600 pairs). Next, for each pair, we varied three additional parameters in the simulations: the spatial extent of the sources, the strength of the coupling between the two sources, and the signal-to-noise ratio (SNR) at the sensors. These parameters are described in more detail below.

*(i) Patches and Point-Like Sources*. We simulated point-like sources (i.e., 1 dipole) and cortical patches (with surface areas: 2, 4, or 8 cm^{2}). The activity of a cortical patch was simulated by placing identical time series at the vertices that make up the patch. As described above, the vertices were obtained via tessellation of the MNI brain (Colin 27).

*(ii) Coupling Strength*. When generating the time series for each of the two sources of a pair of simulated generators, we defined the coupling by setting the alpha-band coherence to either 0.1, 0.2, or 0.4. We simulated time series that were 7000 samples long (600 Hz sampling frequency). Note here that by actually simulating true coherence at the source level, we circumvent debates about spurious coupling that arises from field spread. As such, we assess here the ability to recover truly coherent cortical activity.

*(iii) Signal*-*to*-*Noise Ratio (SNR)*. White noise was added to the sensor signals in order to achieve three levels of SNR (0 dB, −20 dB, and −40 dB), calculated as the ratio of the Frobenius norm of signal and noise amplitudes at the sensors.

In summary, we randomly chose 600 different pairs of cortical location configurations, for which we varied source size (4 levels: point-like, 2, 4, or 8 cm^{2}), coupling strength (3 levels: 0.1, 0.2, or 0.4), and SNR (3 levels: 0 dB, −20 dB, and −40 dB). This yielded a total of 21,600 sets of simulated MEG sensor-level data. Next, we evaluated the effect of varying the Tikhonov regularization parameter on detection performance of MNE, independently, for power and for coherence mapping.

##### 2.4. Power and Coherence Reconstructions

The previous section described the cortical power and coherence configurations that were generated in the MEG data simulation step. Now, in order to assess our ability to reconstruct these simulated (i.e., known) source configurations, we first reconstruct the source time series by applying MNE to the simulated sensor data, and then we compute spectral power and coherence from these estimated time series.

Most MEG studies that use MNE select a single lambda value once and for all to be used throughout the study; the same regularization coefficient is hence used for both power and coupling estimations (assuming both are performed within the study). This is not the case here, precisely because our objective is to test whether optimal lambda values differ between power and coherence mapping.

Therefore, in order to find the regularization coefficient providing the most accurate reconstruction of (a) power and (b) interareal coherence across all 21,600 simulated configurations, we used multiple values for ranging from to (7 values). In addition to searching for the optimal lambda that provides the best results across all simulations, the definition of best lambda was also separately examined as a function of SNR, coupling strength, and source extent. Finally, for the sake of comparison, we also computed the optimal value obtained from the standard L-curve method [51]. The range of possible lambda values examined in our estimations was selected based on visual inspection in a subset of simulations and chosen to ensure it included the lambda obtained via the L-curve approach.

Notably, for the coupling analysis, we used the reconstructed time series at one of the two simulated dipoles as reference point, and we calculated coherence with all other source time series across the brain. Spectral power (at any location ) and magnitude squared coherence (between two locations and ) were calculated as follows:where is the cross-spectral density matrix and is the frequency bin. Power and MS coherence were calculated via built-in standard MATLAB (Mathworks Inc., MA, USA) functions based on Welch’s averaged and modified periodogram method [52].

##### 2.5. Receiver Operating Characteristic (ROC) Curves

We evaluated the performance of each method using the area under the curve (AUC) from the ROC curves obtained by plotting the True Positive Fraction (TPF) versus the False Positive Fraction (FPF) at a given threshold , calculated as follows: where TP() represents the true positives (defined as the intersection between simulated sources and active sources at threshold ) and FP() represents the false positives (defined as all active sources but excluding the true positives at activation threshold ). By computing TPF and FPF repeatedly for successive values of activation thresholds , we obtain a ROC curve from which we derive the AUC. The AUC is taken as a measure of performance; that is, the best regularization coefficient would be the one that yields the highest AUC. Statistical comparisons were performed using standard parametric two-tailed -tests. Note that, for reference-based coherence reconstruction, computing true positives can be ambiguous if we consider the sources within the “reference patch” (location 1) to be true positives. To avoid this problem, we chose to quantify how well the distant coherent patch (location 2) was detected. In other words, we used the time series estimated at location 1 (as a seed) and considered only the simulated activity that make up the patch at location 2 to be the vertices we wish to detect. Therefore, the vertices within the reference patch were excluded from the ROC calculations for coherence evaluations.

#### 3. Results

##### 3.1. Optimal Tikhonov Regularization Coefficient Is Not the Same for Power and Coherence Reconstructions

Figure 1(a) shows the effect of lambda selection on quality of power and coherence detection across all simulated conditions and configurations (measured as mean AUC). Importantly, we found that the best mean value of lambda for power () differs from the best mean lambda for coherence, which was two orders of magnitude smaller (). In comparison to the optimal value for power, the lower optimal value for coherence implies that a smaller residual should be allowed to have an optimal coherence reconstruction. The observations in Figure 1(a) also indicate that coherence reconstructions are more sensitive to the selection of an appropriate lambda value than power reconstructions (as AUC for coherence peaks at and drops again, while AUC for power displays a flatter distribution for values neighboring ).