Minimizing the distortions in electrophysiological source imaging of cortical oscillatory activity via Spectral Structured Sparse Bayesian Learning

Oscillatory processes at all spatial scales and on all frequencies underpin brain function. Electrophysiological Source Imaging (ESI) is the data-driven brain imaging modality that provides the inverse solutions to the source processes of the EEG, MEG, or ECoG data. This study aimed to carry out an ESI of the source cross-spectrum while controlling common distortions of the estimates. As with all ESI-related problems under realistic settings, the main obstacle we faced is a severely ill-conditioned and high-dimensional inverse problem. Therefore, we opted for Bayesian inverse solutions that posited a priori probabilities on the source process. Indeed, rigorously specifying both the likelihoods and a priori probabilities of the problem leads to the proper Bayesian inverse problem of cross-spectral matrices. These inverse solutions are our formal definition for cross-spectral ESI (cESI), which requires a priori of the source cross-spectrum to counter the severe ill-condition and high-dimensionality of matrices. However, inverse solutions for this problem were NP-hard to tackle or approximated within iterations with bad-conditioned matrices in the standard ESI setup. We introduce cESI with a joint a priori probability upon the source cross-spectrum to avoid these problems. cESI inverse solutions are low-dimensional ones for the set of random vector instances and not random matrices. We achieved cESI inverse solutions through the variational approximations via our Spectral Structured Sparse Bayesian Learning (ssSBL) algorithm https://github.com/CCC-members/Spectral-Structured-Sparse-Bayesian-Learning. We compared low-density EEG (10–20 system) ssSBL inverse solutions with reference cESIs for two experiments: (a) high-density MEG that were used to simulate EEG and (b) high-density macaque ECoG that were recorded simultaneously with EEG. The ssSBL resulted in two orders of magnitude with less distortion than the state-of-the-art ESI methods. Our cESI toolbox, including the ssSBL method, is available at https://github.com/CCC-members/BC-VARETA_Toolbox.

Introduction

Brain function is embodied in large-scale dynamic networks underlying all behavior and cognition. The natural modes of these networks are oscillations at functionally specific frequencies (Engel et al., 2001; Varela et al., 2001). Direct (invasive) observations of brain oscillations are increasingly more fine-grained and informative (Buzsáki et al., 2013; Frauscher et al., 2018). However, they cannot be feasibly examined in the whole brain in vivo human studies. Hence, considerable effort is dedicated to developing indirect (non-invasive) imaging methods to reveal the brain's oscillatory architecture with fine-grained time resolution. Notably, such an imaging modality should reveal the two different but interrelated aspects of networks: (a) the activity of each dynamical unit (node) and (b) their connectivity.

Electrophysiological Source Imaging (ESI) is a prime candidate for mapping brain network activity and connectivity. ESI is predicated upon the fact that the Electro-encephalogram (EEG), Magneto-encephalogram (MEG), or Electrocorticogram (ECoG) are generated by the electrical currents of macroscopic neural masses (nodes) resulting from the local mean field of post-synaptic potential activity (Jirsa and Haken, 1996, 1997; Jirsa, 2004; Deco et al., 2008; Friston, 2009; Daunizeau et al., 2010, 2011; Moran et al., 2013).

In addition, ESI has been used in attempts to estimate these currents, also known as source activity, from their electrophysiological observables. For a neural mass, the current is directly proportional to the mean field activity (Valdes-Sosa et al., 2009; Rosa et al., 2010). In turn, each neural mass is the collection of neurons to the extent of a few millimeters such that a mean-field observable (Freeman, 1975; Vinck and Perrenoud, 2019) is the “source” of the observations. In this study, we attempted the most general formulation of the problem as a reference for future work while providing concrete examples.

In principle, source activity in a network is modeled as a strictly dynamic random field ι(ϱ , ς) over a continuous spatiotemporal manifold with ϱ ∈ ℝ³ (the parts of the brain that generate observations) and continuous time ς ∈ ℝ. In practice, the random field ι(ϱ , ς) must be discretized at spatial points ϱ_g; g = 1, ⋯ , G, and at discrete time points ς_t = t △ ς; t = −T, ⋯ , T. In this scenario △ς is the sampling period. Thus, we focused on the vector time series ι(t), defined as the vector function with entries ι(g , t) = ι(ϱ_g , ς_t). The (multi-channel) electro-physiological data is the vector time series v(t) with entries v(e , t) for each sensor, e = 1⋯E that arises from the discretization of v(ϱ , ς), where the electromagnetic field was produced by ι(ϱ, ς).

The data v(t) from its latent source activity ι(t) is presented in the following forward model (Eq. 1).

$\begin{array}{ll}v\left(t\right)={\text{L}}_{v\iota }\iota \left(t\right)+\xi \left(t\right)& (1)\end{array}$

Where L_vι is the real-valued lead field matrix or forward operator that projects sources ι(t) to forward the data v(t), and ξ(t) is the time series of instrumental noise assumed to be independent of the source activity ι(t). The forward operator (lead field) L_vι is linear and stationary by definition, derived from the discretization of a quasi-static electromagnetic forward model (Hämäläinen et al., 1993; Riera and Fuentes, 1998; Hallez et al., 2007; Lei et al., 2011; Piastra et al., 2020). For operators, we used suffixes that indicate the operator's codomain and domain.

ESI can be defined as the generally non-linear inverse solution (Eq. 2) (Nunez, 1974; Hämäläinen and Ilmoniemi, 1994; Nunez et al., 1994; Baillet et al., 2001; Nunez and Srinivasan, 2006; Burle et al., 2015) via the optimal inverse operator that we denote with a hat ${\widehat{\text{W}}}_{\iota v}$. An inverse operator ${\widehat{\text{W}}}_{\iota v}$ projects the data v(t) to explain approximately its source and produce its estimator $\widehat{\iota }\left(t\right)$.

$\begin{array}{ll}\widehat{\iota }\left(t\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left(v\left(t\right)\right)& (2)\end{array}$

Optimizing W_ιv from the data solves an inverse problem, Eq. (1), that is not only ill-posed in the sense of Hadamard (Hadamard and Morse, 1953) with degeneracy in a (G − S)-dimensional space but also severely ill-conditioned and high dimensional (with G ≫ E). Overcoming these challenges to obtain acceptable inverse solutions has been the subject of much research in specific optimization methods for W_ιv which was well-summarized in the study of Knösche and Haueisen (2022).

Although we have defined ESI in terms of time domain signals, it is well-known that brain activity is oscillatory at all scales, from the local field potentials at the neuronal level (Freeman, 1975; Vinck and Perrenoud, 2019) to the EEG, MEG, or ECoG (Niedermeyer and da Silva, 2005; Le Van Quyen and Bragin, 2007; Frauscher et al., 2018). In tailoring ESI for oscillatory activity, a natural framework is that of the frequency domain ι(f), a random vector representing the (discrete) Fourier transform of vector time series ι(t), and comprising complex-valued entries ι(g , f) for each source g = 1, ⋯ , G and frequency f = −T, ⋯ , T. Considering that we may compute the physical frequency as ν_f = f△ν for a spectral period $△\nu =\frac{1}{\left(\left(2T+1\right)△\varsigma \right)}$ where (2 T + 1)△ς is the Nyquist frequency, the corresponding frequency domain data term is v(f).

Then, the equivalent expressions to the previous Eqs. (1) and (2) in the frequency domain are Eq. (3), the corresponding inverse problem for the Fourier transform ι(f), and its inverse solution leading to the estimator $\widehat{\iota }\left(f\right)$. Considering that solving an inverse problem in the frequency domain should be (ideally) optimizing the inverse operator ${\widehat{\text{W}}}_{\iota v}$ Eq. (3) from the data Fourier transform v(f) (and not from the process v(t)).

$\begin{array}{ll}\begin{array}{c}v\left(f\right)={\text{L}}_{v\iota }\iota \left(f\right)+\xi \left(f\right)\\ \widehat{\iota }\left(f\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left(v\left(f\right)\right)\end{array}& (3)\end{array}$

Frequency domain descriptions of the electrophysiological observations and their inverse solutions (Eqs. 1–3) have been used fruitfully to probe behavior and cognition (Valdés et al., 1992; Engel et al., 2001; Varela et al., 2001; Le Van Quyen and Bragin, 2007; Marzetti et al., 2008; Valdes-Sosa et al., 2009; Brookes et al., 2011a,b; Faes and Nollo, 2011; Faes et al., 2012, 2017; Friston et al., 2012; Hipp et al., 2012; Colclough et al., 2015; Wens et al., 2015; Mahjoory et al., 2017; Vidaurre et al., 2018a,b; Tewarie et al., 2019; Nolte et al., 2020).

Our concern, then, is with frequency domain ESI in Eq. (3), our primary target being the source cross-spectral density matrix or source cross-spectrum ${\Sigma }_{\iota \iota }\left(f\right)=\langle \iota \left(f\right){\iota }^{†}\left(f\right)\rangle$, with the expected value over the sample space of the discrete Fourier transform ι(f) (Valdés et al., 1992; Engel et al., 2001; Varela et al., 2001; Nunez and Srinivasan, 2006; Hipp et al., 2012; Vidaurre et al., 2018b). The corresponding data cross-spectrum is ${\Sigma }_{vv}\left(f\right)=\langle v\left(f\right)v^{†}\left(f\right)\rangle$. In practice, this later quantity must be substituted by its estimator ${\overline{\Sigma }}_{vv}\left(f\right)$ Eq. (4).

$\begin{array}{ll}{\overline{\Sigma }}_{vv}\left(f\right)=\langle v_m\left(f\right)v_m^{†}\left(f\right);M\rangle =\frac{1}{M}{\displaystyle \sum _{m=1}^M}{v}_m\left(f\right)v_m^{†}\left(f\right)& (4)\end{array}$

${\overline{\Sigma }}_{vv}\left(f\right)$ is denoted with a different hat type since the expectation is for a finite number of instances M with an index m = 1⋯M. Toward this estimation, we followed the standard practice of segmenting the data time series v(t) into segments (v_m (t); ∀ m). Thus, we worked with instances (v_m (f); ∀ m) of the discrete Fourier transform applied to realizations or observations for these segments.

The frequency domain inverse problem and the inverse solution for the source cross-spectrum Σ_ιι(f) may be stated as Eq. (5) (He et al., 2019), valid under the condition of independence between the source process ι(t) and the noise process ξ(t) in the previous forward model (Eq. 1).

$\begin{array}{ll}\begin{array}{c}{\Sigma }_{vv}\left(f\right)={\text{L}}_{v\iota }{\Sigma }_{\iota \iota }\left(f\right){\text{L}}_{\iota v}+{\Sigma }_{\xi \xi }\left(f\right)\\ {\widehat{\Sigma }}_{\iota \iota }\left(f\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left({\overline{\Sigma }}_{vv}\left(f\right)\right)\\ \end{array}& (5)\end{array}$

The pursuit of the cross-spectrum Σ_ιι(f) is rewarding since it completely specifies the multivariate linear properties of any stochastically driven system, be it linear or non-linear (Brillinger, 2001, 2012), though the latter requires additional higher-order kernels for a complete description (Brillinger, 1965; Brillinger and Rosenblatt, 1967). We used the asymptotic stochastic properties of the discrete Fourier transform to introduce these developments (Section Asymptotic probability theory of the Fourier transform and cross-spectrum).

Cross-spectral diagonal elements Σ_ιι(g , g , f) have intuitive interpretations: the variances ${\sigma }_{\iota \iota }^2\left(g,f\right)$ (${\sigma }_{\iota \iota }^2\left(f\right)=diag\left({\Sigma }_{\iota \iota }\left(f\right)\right)$) are the spectra of source activity, the Cortical Spectral Topography (CST). Off-diagonal elements ${\Sigma }_{\iota \iota }\left(g,g{}^{\prime },f\right)$ are the cross-spectra reflecting functional connectivity. Optimal inverse operators ${\widehat{\text{W}}}_{\iota v}$ for the cross-spectrum Σ_ιι(f) Eq. (5) is a novel form of electrophysiological source imaging: cross-spectral ESI (cESI).

From expressions (Eqs. 1–5), the time and frequency domain variants of the inverse problem, it is evident that cESI, the inverse solution ${\widehat{\Sigma }}_{\iota \iota }\left(f\right)$ of the source cross-spectrum Σ_ιι(f), may be obtained from three different types of primary information: v(t), or v(f), or ${\overline{\Sigma }}_{vv}\left(f\right)$. Initially, it might seem that an inverse solution ${\widehat{\Sigma }}_{\iota \iota }\left(f\right)$ from any of these data types would be equivalent. However, as shown further into the study, each data type requires optimizing its specific inverse operator ${\widehat{\text{W}}}_{\iota v}$, thus defining different estimation “routes” to cross-spectrum Σ_ιι(f).

We describe in detail the theory, benefits, and problems of each route, in terms of the forward operator L_vι (Section Forward projection by the lead field linear and stationary operator), and the inverse operator W_ιv (Section Bayesian (MAP) inverse operators. Quasilinear (F-invariant) approximation).

In order to explore these routes in the following sections (Sections Bayesian MAP1 inverse operators. MNE, eLORETA, and LCMV as particular cases and Data used for the validation), we need to select an inverse solution framework, of which there are many potentially valuable approaches (Knösche and Haueisen, 2022). We adopted the Bayesian maximum a posteriori (MAP) probability approach (MacKay, 2003). The MAPs (for each route) are derived from finding the latent variables ${X}$ that maximize the posterior probability $q\left({X}|{Y}\right)$(Eq. 6).

$\begin{array}{ll}q\left({X}|{Y}\right)\propto p\left({Y}|{X}\right)p\left({X}\right)& (6)\end{array}$

Where $p\left({Y}|{X}\right)$ is the likelihood of the data ${Y}$ conditional on the latent variable ${X}$. In our context, ${Y}$ can corresponding to be v(t), or v(f), or ${\overline{\Sigma }}_{vv}\left(f\right)$, and ${X}$ can be corresponding to ι(t), or ι(f), or Σ_ιι(f) (Grave de Peralta Menendez et al., 2004; Friston et al., 2006; Mattout et al., 2006; Friston K. J. et al., 2008; Wipf and Nagarajan, 2009). The term $p\left({X}\right)$ is the respective a priori probability upon the latent variable ι(t), or ι(f), or Σ_ιι(f).

As can be observed in the following section (Section Data used for the validation), the third route inverse solution (Eq. 5) has desirable properties for cESI. However, if not impossible in a realistic cESI setup, such inverse solutions are N-P hard and therefore require Approximated Bayesian Computation (ABC) (Csilléry et al., 2010). This issue has been discussed in detail in the Bayesian literature (Dempster et al., 1977; Liu and Rubin, 1994; Daunizeau and Friston, 2007; Friston et al., 2007; Nummenmaa et al., 2007; Friston K. J. et al., 2008; Paz-Linares et al., 2018).

In this study, we took a more practical path. Rather than attempting to use complicated ABCs to obtain more general solutions, the cESI rationale we focused on is to restrict Bayesian inverse operators W_ιv in Eqs. (1–5) within the “quasilinear” class W_ιv ≅ T_ιv. Quasilinear inverse operators are also known in a more general context as the linear proximal operators or as linear back projectors that solve non-linear optimization or inverse problems (Kaplan and Tichatschke, 1998; Piotrowski and Yamada, 2008; Gramfort et al., 2012; Tirer and Giryes, 2020).

Here, quasilinear inverse operator T_vι, which holds the same linearity attribute as the forward operator L_vι, produced cESI that preserved the amplitude and phase information in the frequency domain (Mantini et al., 2007; Marzetti et al., 2008; Brookes et al., 2011b; Hipp et al., 2012; Lopes da Silva, 2013; He et al., 2019; Nolte et al., 2020). That is, we ensured that T_vι possesses the attribute of what we call “F-invariance” for essential properties (Brillinger, 2001, 2012).

An immediate consequence of taking W_ιv ≅ T_ιv, a quasilinear approximation of the inverse operator in the third cESI route (Eq. 5), is the following representation of source cross-spectrum ${\Sigma }_{\iota \iota }\left(f\right)\cong {\overline{\Sigma }}_{\iota \iota }\left(f\right)$ (Eq. 7). Here, we considered ${\overline{\Sigma }}_{\iota \iota }\left(f\right)$ to be calculated from the set of random instances (ι_m (f); ∀ m) of the latent vector process ι(f).

$\begin{array}{ll}{\overline{\Sigma }}_{\iota \iota }\left(f\right)=\langle {\iota }_m\left(f\right){\iota }_m^{†}\left(f\right);M\rangle =\frac{1}{M}{\displaystyle \sum _{m=1}^M}{\iota }_m\left(f\right){\iota }_m^{†}\left(f\right)& (7)\end{array}$

Once ${\Sigma }_{\iota \iota }\left(f\right)\cong {\overline{\Sigma }}_{\iota \iota }\left(f\right)$ was assumed, the ABCs were simplified dramatically since a MAP for the matrix ${\overline{\Sigma }}_{\iota \iota }\left(f\right)$ (Eq. 7) turns into a joint-MAP (Section Validation rationale) for the random vectors (ι_m (f); ∀ m) that are implicit in ${\overline{\Sigma }}_{\iota \iota }\left(f\right)$ (Hsiao et al., 1998, 2002; Yeredor, 2000; Davis et al., 2001; Auranen et al., 2005; Wipf and Nagarajan, 2009; Chen et al., 2011).

We implemenedt the joint-MAP via Spectral Structured Sparse Bayesian Learning (ssSBL), the type of ABC developed in Section Measures of distortion. As shown under realistic inverse problem settings (Section Results), the third cESI route (Eq. 5) implemented via the ssSBL approach leads to less distorted estimates than the traditional methods for the first and second cESI routes (Eqs. 1–3). To judge distortions, we employed the well-known ESI methods as a baseline: Exact Low-Resolution Electromagnetic Tomographic Analysis (ELORETA) (Pascual-Marqui et al., 2006) and Linearly Constrained Minimum Variance (LCMV) (Van Veen et al., 1997).

The theoretical framework allowed one to consider the fundamental problem of ESI distortions. Indeed, significant distortions are expected with any state-of-the-art inverse solutions in a realistic ESI setup. The distortions, which we explore later, are localization error and leakage (blurring). These distortions are pervasive comparing simulated topographic vectors, say ι(t) or ι(f), vs. their inverse solution $\widehat{\iota }\left(t\right)$ or $\widehat{\iota }\left(f\right)$ (Kobayashi et al., 2003; Grova et al., 2008; Schoffelen and Gross, 2009; Haufe et al., 2013; Burle et al., 2015; Colclough et al., 2015; Bradley et al., 2016; Mahjoory et al., 2017; Stokes and Purdon, 2017; He et al., 2018, 2019; Palva et al., 2018; Haufe and Ewald, 2019; Marinazzo et al., 2019).

As we will show in Section Discussion, the topographic distortions of the inverse solutions for a random vector ι(t) or ι(f) can reach unacceptable levels for second-order statistics, such as the sample estimator for the cross-spectrum ${\overline{\Sigma }}_{\widehat{\iota }\widehat{\iota }}\left(f\right)$ calculated from an inverse solution $\widehat{\iota }\left(t\right)$ or $\widehat{\iota }\left(f\right)$. Minimizing CST distortions (for the estimator of the spectrum ${\sigma }_{\iota \iota }^2\left(f\right)$) benefits the overall cross-spectral estimation (for Σ_ιι(f)). Our results suggest that that ssSBL significantly reduces these distortions.

Standard cESI theory

Asymptotic probability theory of the Fourier transform and cross-spectrum

The material in this section (with minor differences in notation) is described in greater detail in Brillinger (2012). For an exhaustive definition of terms, refer to Supplementary material (SD, Section Introduction). Our primary interest will be in frequency domain quantities. Let x(t) be an R-dimensional vector time series. We worked with the following definitions:

• The vector stochastic process x(f) (Eq. 8) is defined in the discrete frequency domain ν_f = f △ ν with f = −T, ⋯ , T, as the discrete Fourier transform of the vector time series x(t).

$\begin{array}{ll}\text{x}\left(f\right)={\displaystyle \sum _{t=-T}^T}\text{x}\left(t\right)e^{-i2\pi \left(f△\nu \right)\left(t△\varsigma \right)}& (8)\end{array}$

• The cross-spectral density matrix or cross-spectrum Σ_xx(ν) (Eq. 9) was defined in the frequency domain ν as the Fourier transform of the auto-covariance matrix Σ_xx(τ).

$\begin{array}{ll}{\Sigma }_{\text{x}\text{x}}\left(\nu \right)={\displaystyle \sum _{\tau =-\infty }^{+\infty }}{\Sigma }_{\text{x}\text{x}}\left(\tau \right)e^{-i2\pi \nu \left(\tau △\varsigma \right)}& (9)\end{array}$

Here, the auto-covariance ${\Sigma }_{\text{x}\text{x}}\left(\tau \right)=\langle \text{x}\left(t\right){\text{x}}^{†}\left(t+\tau \right)\rangle$ depends on the time-lag τ and does not vary with time t, thus being second-order stationary. We will furthermore assume (for technical reasons) that x(t) is a strictly stationary vector time series where all moments are also translation invariants. We also assumed that the strong mixing condition holds. This condition was due to the rapid decrease in the magnitude of the autocovariance Σ_xx(τ), and all higher-order moments as the time-lag τ increases.

A fundamental result on which we based our work is Theorem 4.1.1 of Brillinger (2001), which can be understood as the equivalent for Fourier coefficients of the central limit theorem under stationarity and strong mixing conditions (Rosenblatt, 1956). In our notation, this theorem states:

“Assume that the number of time points in the discrete Fourier transform (Eq. 8), goes to the infinity (T → +∞), and let the sampling period go to zero (△ς → 0+) so that the spectral resolution $△\nu =\frac{1}{\left(\left(2T+1\right)△\varsigma \right)}=constant$ holds constant. Then, it holds that x(f) is asymptotically independent for all f and converges in probability to the circularly symmetric multivariate complex-valued Gaussian probability density (Eq. 10). Where the Hermitian covariance matrix Σ_xx(f) is the cross-spectrum (Eq. 9) at the frequencies ν_f = f △ ν with f = −T, ⋯ , T.”

$\begin{array}{ll}{N}^{ℂ}\left(\text{x}\left(f\right)|0,{\Sigma }_{\text{x}\text{x}}\left(f\right)\right)=\frac{1}{\left|\pi {\Sigma }_{\text{x}\text{x}}\left(f\right)\right|}{e}^{-{\text{x}}^{†}\left(f\right){\Sigma }_{\text{x}\text{x}}^{-1}\left(f\right)\text{x}\left(f\right)}& (10)\end{array}$

We emphasize that this Gaussian distribution asymptotic distribution not only holds for x(f) but also for time-varying estimators such as the time-windowed discrete Fourier transform or the Hilbert transform x(t , f) (Bruns, 2004). The latter has been widely used in the literature (Faes and Nollo, 2011; Friston et al., 2012; Nolte et al., 2020).

It is important to note that Gaussianity might not be valid for the original time series x(t) or even one of its band-filtered versions, as many assume in the literature (Grave de Peralta Menendez et al., 2004; Friston et al., 2006; Mattout et al., 2006; Friston K. J. et al., 2008; Wipf and Nagarajan, 2009; Paz-Linares et al., 2017). In fact, even in the case of non-Gaussian, non-linear time series Brillinger's theorem is valid as long as stationarity and strong mixing hold. This validity does not imply that the spectral density matrix Σ_xx(f) completely characterizes the non-linear or non-Gaussian system. Cumulant information of orders higher than two may be necessary for a complete system description (Brillinger, 1965; Brillinger and Rosenblatt, 1967).

Brillinger's theorem leads to the probability density of any sampled estimator of the cross-spectrum. In particular, it applies to ${\overline{\Sigma }}_{\text{x}\text{x}}\left(f\right)=\langle {\text{x}}_m\left(f\right){\text{x}}_m^{†}\left(f\right);M\rangle$ the sampled estimator for Fourier transform instances (x_m (f); ∀ m) with sample size M. Here, these instances (x_m (f); ∀ m) represent the discrete Fourier transform applied to sample realizations (x_m (t); ∀ m) obtained from time segments of the observations. Then, it follows that a Hermitian Wishart W^ℂ (Eq. 11), with R-dimensional scale matrix Σ_xx(f) (cross-spectrum), and degree of freedom M (with M ⩾ R), is the probability density of the estimator ${\overline{\Sigma }}_{\text{x}\text{x}}\left(f\right)$.

$\begin{array}{ll}{W}^{ℂ}\left({\overline{\Sigma }}_{\text{x}\text{x}}\left(f\right)|{\Sigma }_{\text{x}\text{x}}\left(f\right),M\right)\propto \frac{|{\overline{\Sigma }}_{\text{x}\text{x}}\left(f\right){|}^{M-R}}{|{\Sigma }_{\text{x}\text{x}}\left(f\right){|}^M}{e}^{-Mtr\left({\Sigma }_{\text{x}\text{x}}^{-1}\left(f\right){\overline{\Sigma }}_{\text{x}\text{x}}\left(f\right)\right)}& (11)\end{array}$

Since we based our further developments on the assumption of Gaussianity (Eqs. 10, 11), we carried out the statistical test for the distribution of the discrete Fourier transform x(f) (Eq. 10) with two examples of resting state sensor data: MEG data from the Human Connectome Project (HCP) (Van Essen et al., 2013) and Macaque ECoG data from the Neurotycho project (Nagasaka et al., 2011). The outcome of this test for data did not allow us to reject the hypothesis of Gaussianity, which then was also plausible for the sources.

Forward projection by the lead field linear and stationary operator

We emphasize that linearity and stationarity are essential “conservative” attributes of the operators for our target (cESI). We note, this is the reason why one can state all the forward “routes” in terms of the same operator L_vι (Eq. 12). Preserving the Gaussianity of the Fourier transform ι(f) is only possible under the linear forward operator L_vι and subsequent linear or quasilinear inverse operator T_ιv (Marzetti et al., 2008; He et al., 2019; Nolte et al., 2020).

$\begin{array}{ll}\begin{array}{c}\left(route1\right)\text{v}\left(t\right)={\text{L}}_{v\iota }\iota \left(t\right)+\xi \left(t\right)\\ \left(route2\right)\text{v}\left(f\right)={\text{L}}_{v\iota }\iota \left(f\right)+\xi \left(f\right)\\ \left(route3\right){\Sigma }_{vv}\left(f\right)={\text{L}}_{v\iota }{\Sigma }_{\iota \iota }\left(f\right){\text{L}}_{\iota v}+{\Sigma }_{\xi \xi }\left(f\right)\\ \end{array}& (12)\end{array}$

Furthermore, both attributes are crucial to avoid non-linear warping and delays of the amplitude and phase information in the frequency domain for cESI (Reid et al., 2019). Therefore, we define operators with these properties for the frequency domain as “F-invariant.” Though we restricted our attention later to F-invariant operators T_ιv, non-linear operators W_ιv have been useful in other contexts (Picton and Hillyard, 1974; Picton et al., 1974; Lopes da Silva et al., 1991; Clark et al., 1994; Makeig et al., 1999, 2004; Makeig, 2002; Eichele et al., 2005; Harrison et al., 2008; Vega-Hernández et al., 2008; Maurer and Dierks, 2012).

Noteworthily, distinguishing forward “routes” leads to variants of the inverse problems or inverse operators W_ιv for estimating the specific latent variable ι(t), or ι(f), or Σ_ιι(f), as discussed in the next section (Section Bayesian (MAP) inverse operators. Quasilinear (F-invariant) approximation). Estimation of the cross-spectrum Σ_ιι(f), involves the challenging inverse problem of the matrix equation (route 3) detailed in section (Section Bayesian MAP1 inverse operators. MNE, eLORETA, and LCMV as particular cases).

We currently illustrate the effects of the forward operator with the topographic maps: Cortical Spectral Topography (CST), a map of the cortical spectrum ${\overline{\sigma }}_{\iota \iota }^2\left(f\right)=diag\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)$, and Sensor Spectral Topography (SST), a map of the sensor spectrum ${\overline{\sigma }}_{vv}^2\left(f\right)=diag\left({\overline{\Sigma }}_{vv}\left(f\right)\right)$ (Figure 1). On the left is a hypothetical CST, and on the right is the corresponding SST. The inverse problem consists in estimating the latent CST from observed SST.

FIGURE 1

Figure 1. An illustration of Σ_ιι(f) the source cross-spectrum of human cortical alpha activity and its EEG cross-spectrum Σ_vv(f). For this illustrations we calculate the sampled estimates ${\overline{\Sigma }}_{\iota \iota }\left(f\right)$ and ${\overline{\Sigma }}_{vv}\left(f\right)$ of which we represent the topographic projections: (a) Cortical Spectral Topography (CST) ${\overline{\sigma }}_{\text{u}}^2\left(f\right)={\overline{\Sigma }}_{\text{u}}\left(f\right)$ and (b) Sensor Spectral Topography (SST) ${\overline{\sigma }}_{vv}^2\left(f\right)={\overline{\Sigma }}_{vv}\left(f\right)$.

Bayesian (MAP) inverse operators, quasilinear (F-invariant) approximation

From the theory of inverse problems (Tarantola, 2005), one may seek inverse solutions for each of the latent variables in the routes corresponding to ι(t), or ι(f), or Σ_ιι(f) (Eq. 13) and Figure 2. In each route, the theoretical inverse operator ${\widehat{\text{W}}}_{\iota v}$ represents the symbolic projection to the source space and these inverse solutions. In turn, inverse operators are determined (analytically or numerically) by data-driven methods that depend on the routes from the data variable: v(t), or v(f), or ${\overline{\Sigma }}_{vv}\left(f\right)$, and the forward operator L_vι.

$\begin{array}{ll}\begin{array}{c}\left(route1\right){\widehat{\text{W}}}_{\iota v}\left(v\left(t\right)\right)\\ \left(route2\right)\widehat{\iota }\left(f\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left(v\left(f\right)\right)\\ \left(route3\right){\widehat{\Sigma }}_{\iota \iota }\left(f\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left({\overline{\Sigma }}_{vv}\left(f\right)\right)\\ \end{array}& (13)\end{array}$

FIGURE 2

Figure 2. Inverse operators from (a1) the MEG/EEG/ECoG sensor signal v_t(t), defined as 3D tensor with every trial and time-point observation, to (b3) the source cross-spectrum at every frequency Σ_ιι(f), via the different cross-spectral Electrophysiological Source Imaging (cESI) routes. Route 1 (a1→b1→b2→b3) via an inverse operator W_ιv in the time-domain first determines (b1) the source processes ι₃(t) and then compute (b2) their Fourier transformι₃(f). Route 2 (a1→a2→b2→b3) first computes (a2) the data Fourier transform v_m(f) and then via an inverse operator W_ιv in the frequency domain determines (b2) their source Fourier transform ι_m(f). Route 3 (a1→a2→a3→b3) first computes (a2) the data Fourier transform v₃(f) and then determines (a3) the data cross-spectrum Σ_vv(f). We revindicate Route 3, which is via an inverse operator W_ιv in the spectral-domain with a priori probabilities upon the source cross-spectrum.

Let any of the data variables be represented generically by ${Y}$, and likewise, ${X}$ represent all types of latent variables in Eq. (13). A maximum a posteriori (MAP) Bayesian inverse operator ${\widehat{\text{W}}}_{{XY}}$ (Eq. 14) achieves the maximum of a posteriori probability $q\left({X}|{Y}\right)$. We define the a posteriori $q\left({X}|{Y}\right)$ as the conditional probability determined from the likelihood $p\left({Y}|{X}\right)$ and a priori $p\left({X}\right)$. Depending on $q\left({X}|{Y}\right)$ an inverse operator ${\widehat{\text{W}}}_{{XY}}$ can be non-linear or linear, computed numerically or analytically, also intractable or tractable.

$\begin{array}{ll}\begin{array}{c}{X}\leftarrow {\widehat{\text{W}}}_{{XY}}\left({Y}\right)\\ {\widehat{\text{W}}}_{{XY}}=argma{x}_{{\text{W}}_{{XY}}}\left(q\left({\text{W}}_{{XY}}\left({Y}\right)|{Y}\right)\right)\\ q\left({X}|{Y}\right)\propto p\left({Y}|{X}\right)p\left({X}\right)\end{array}& (14)\end{array}$

The essential role of the a priori $p\left({X}\right)$ is to overcome the ill condition of the likelihood $p\left({Y}|{X}\right)$, in other words, to provide a unique solution to the inverse problem for the source variable ${X}$. Here we formulate this a priori $p\left({X}\right)$ as a Gibbs probability density (Eq. 15).

$\begin{array}{ll}p\left({X}\right)\propto exp\left(H\left({X}\right)|\alpha \right),& (15)\end{array}$

Where α is the Gibbs temperature parameter. The Gibbs energy function $H\left({X}\right)$ is commonly defined using the norms of vectors or matrices (Petersen and Pedersen, 2008; Golub and Van Loan, 2013), usually used to reflect empirical criteria about the structure and density of the source variable ${X}$ defined over the spatial, time, or frequency domains. Assimilating these source qualities requires, in addition, the empirically determined (data-driven) scale parameter α, by some method from the data variable ${Y}$.

We emphasize that this Bayesian MAP formalism is particularly helpful in developing the cESI routes and our optimal inverse solution in the following sections (Sections Bayesian MAP1 inverse operators. MNE, eLORETA, and LCMV as particular cases and Data used for the validation). In another context, Eqs. (14, 15) may be completely equivalent to the classical Tikhonov regularization (Tikhonov and Arsenin, 1977), taking a posteriori $p\left({X}|{Y}\right)$ under the −log transformation. However, the Bayesian formalism is more general than the Tikhonov regularization since the Gibbs energy $H\left({X}\right)$ (Eq. 15) describes the statistical properties of any physical system (Landau and Lifshitz, 1980).

As stated before, we constrained all routes (Eq. 13) to the class of quasilinear inverse operators T_ιv (Eq. 16) (Baillet et al., 2001; Grech et al., 2008). Moreover, the quasilinear class leads to the type of F-invariant inverse solutions for the properties in Eqs. (10, 11), and particularly the tractable cESI route 3.

$\begin{array}{ll}\begin{array}{c}\left(route1\right)\widehat{\iota }\left(t\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left(v\left(t\right)\right)\cong {\widehat{\text{T}}}_{\iota v}\left(v\left(t\right)\right)v\left(t\right)\\ \left(route2\right)\widehat{\iota }\left(f\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left(v\left(f\right)\right)\cong {\widehat{\text{T}}}_{\iota v}\left(v\left(f\right)\right)v\left(f\right)\\ \left(route3\right){\widehat{\Sigma }}_{\iota \iota }\left(f\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left({\overline{\Sigma }}_{vv}\left(f\right)\right)\cong {\widehat{\text{T}}}_{\iota v}\left({\overline{\Sigma }}_{vv}\left(f\right)\right)\\ {\overline{\Sigma }}_{vv}\left(f\right){\widehat{\text{T}}}_{v\iota }\left({\overline{\Sigma }}_{vv}\left(f\right)\right)\\ \end{array}& (16)\end{array}$

The term “quasilinear” T_ιv is applies to inverse operators defined as non-linear matrix functions of the vector argument v(t) or v(f), or the matrix argument ${\overline{\Sigma }}_{vv}\left(f\right)$. A matrix function comprises entries T_ιv(g , s) in the G × E cartesian product of sources g = 1⋯G and sensors e = 1⋯E. These entries are non-linear functions of the vector arguments T_ιv(g , s , v (t)), or T_ιv(g , s , v (f)), or the matrix argument ${\text{T}}_{\iota v}\left(g,s,{\overline{\Sigma }}_{vv}\left(f\right)\right)$.

Bayesian MAP1 inverse operators. MNE, eLORETA, and LCMV as particular cases

Inverse solutions for routes 1 and 2 (Eq. 13) are the traditional activation Electrophysiological Source Imaging (aESI). We obtained these inverse solutions as a first-type MAP (MAP1) (Eq. 17), which estimates the source variable ι(t) in the time domain (Hämäläinen and Ilmoniemi, 1994), or the source variable ι(f) in the frequency domain (Salmelin and Hämäläinen, 1995).

$\begin{array}{ll}\begin{array}{c}\left(route1\right)\widehat{\iota }\left(t\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left(v\left(t\right)\right)\\ {\widehat{\text{W}}}_{\iota v}=argma{x}_{{\text{W}}_{\iota v}}q\left({\text{W}}_{\iota v}\left(v\left(t\right)\right)|v\left(t\right)\right)\\ q\left(\iota \left(t\right)|v\left(t\right)\right)\propto q\left(v\left(t\right)|\iota \left(t\right)\right)p\left(\iota \left(t\right)\right)\\ \left(route2\right)\widehat{\iota }\left(f\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left(v\left(f\right)\right){\widehat{\text{W}}}_{\iota v}=argma{x}_{{\text{W}}_{\iota v}}q\left({\text{W}}_{\iota v}\left(v\left(f\right)\right)|v\left(f\right)\right)\\ q\left(\iota \left(f\right)|v\left(f\right)\right)\propto p\left(v\left(f\right)|\iota \left(f\right)\right)p\left(\iota \left(f\right)\right)\end{array}& (17)\end{array}$

The a posteriori probabilities require definitions of the likelihood, and the a priori is given below in Eq. (18). For the data v(t) a real-valued Gaussian is commonly assumed with mean L_vιι(t) and noise covariance Σ_ξξ(t). This assumption might not be valid for all types of time domain data. In contrast, and in virtue of Brillinger's theorem cited above, data v(f) obtained with the discrete Fourier transform will almost certainly have a complex-valued Gaussian likelihood with mean L_vιι(f) and noise cross-spectrum Σ_ξξ(f) and thus the one we adopted. The a priori probabilities in Eq. (18) are placed upon the activation Gibbs energy in the time domain H(ι(t)) or the frequency domain H(ι(f)) (Eq. 18) as expressed by the vector p-norm (frequently p-norm^p) upon real-valued ι(t) or complex-valued ι(f) (Riesz, 1910; Rudin, 1970; Dunford and Schwartz, 1988; Bourbaki, 2013).

$\begin{array}{ll}\begin{array}{c}\left(route1\right)p\left(v\left(t\right)|\iota \left(t\right)\right)=N^{ℝ}\left(v\left(t\right)|{\text{L}}_{v\iota }\iota \left(t\right),{\Sigma }_{\xi \xi }\left(t\right)\right)\\ p\left(\iota \left(t\right)\right)\propto exp\left(H\left(\iota \left(t\right)\right)|\alpha \left(t\right)\right)\\ \left(route2\right)p\left(v\left(f\right)|\iota \left(f\right)\right)=N^{ℂ}\left(v\left(f\right)|{\text{L}}_{v\iota }\iota \left(f\right),{\Sigma }_{\xi \xi }\left(f\right)\right)\\ p\left(\iota \left(f\right)\right)\propto exp\left(H\left(\iota \left(f\right)\right)|\alpha \left(f\right)\right)\\ \end{array}& (18)\end{array}$

Henceforth, we ruled out route 1 since it is based on the ad-hoc Gaussian assumption for p(v(t)|ι(t)), which, as mentioned before, is not always tenable. Furthermore, our interest was in the frequency domain, which concentrates our attention on route 2, and which bases the likelihood of the Fourier transform p(v(f)|ι(f)) (Eq. 18) on the complex-valued Gaussian probability (Eq. 10).

Toward cESI, selecting a priori p(ι(f)) (Eq. 18) follows the standard rationale that oscillatory brain networks and their activity are characterized by a large-scale and dense (or non-sparse) distribution in conditions of resting-state or task in a block design (Mantini et al., 2007; Brookes et al., 2011b; Hipp et al., 2012; Lopes da Silva, 2013). Such activity is the stochastic and stationary process ι(t) that is composited by oscillations in the frequency domain, as described by the Fourier transform ι(f) (Engel et al., 2001; Varela et al., 2001; Vidaurre et al., 2018a,b; Tewarie et al., 2019).

A general smooth a priori model posits the following Gibbs energy H_{2, A(f)}(ι(f)) (Eq. 19) which, in addition, specifies A(f) as some positive definite and symmetric, or Hermitian, matrix. Specifying the matrix A(f) may follow some types of goodness criteria of the inverse solution, and data-driven methods, which lead to the most common cases of quasilinear (F-invariant) inverse operators (Baillet et al., 2001; Hauk, 2004; Friston K. J. et al., 2008; Grech et al., 2008; Marzetti et al., 2008; Henson et al., 2011; Hindriks, 2020).

$\begin{array}{ll}{H}_{2,\text{A}\left(f\right)}\left(\iota \left(f\right)\right)={\iota }^{†}\left(f\right){\text{A}}^{-1}\left(f\right)\iota \left(f\right)& (19)\end{array}$

Representing the most common quasilinear cases the following inverse operator ${\widehat{\text{T}}}_{\iota v}\left(\text{A}\left(f\right),\text{B}\left(f\right)\right)$, in brief notation ${\widehat{\text{T}}}_{\iota v}\left(f\right)$ (Eq. 20), the constrained generalized inverse of the forward operator L_ιv that incorporates the regularization matrices A(f) and B(f).

$\begin{array}{ll}\begin{array}{c}\widehat{\iota }\left(f\right)\leftarrow {\widehat{\text{T}}}_{\iota v}\left(f\right)v\left(f\right)\\ {\widehat{\text{T}}}_{\iota v}\left(f\right)={\widehat{\Pi }}_{\iota \iota }\left(f\right){\text{L}}_{\iota v}{\text{B}}^{-1}\left(f\right)with{\widehat{\Pi }}_{\iota \iota }\left(f\right)=({\text{L}}_{\iota v}{\text{B}}^{-1}\left(f\right){\text{L}}_{v\iota }\\ +\alpha \left(f\right){\text{A}}^{-1}\left(f\right){)}^{-1}\\ q\left(\iota \left(f\right)|v\left(f\right)\right)\propto N^{ℂ}\left(\iota \left(f\right)|{\widehat{\text{T}}}_{\iota v}\left(f\right)v\left(f\right),{\widehat{\Pi }}_{\iota \iota }\left(f\right)\right)\\ \end{array}& (20)\end{array}$

Where ${\widehat{\text{T}}}_{\iota v}\left(f\right)$ arises from the route 2 MAP1 (Eq. 17) that incorporates the Gibbs energy H_{2, A(f)}(ι(f)) (Eq. 19), and B(f) any approximation to the noise cross-spectrum Σ_ξξ(f) in probabilities (Eq. 18); and q(ι(f)|v(f)) is the Gaussian a posteriori N^ℂ, with mean ${\widehat{\text{T}}}_{\iota v}\left(f\right)v\left(f\right)$ and covariance matrix ${\widehat{\Pi }}_{\iota \iota }\left(f\right)$, that follows from the conjugated relation between the likelihood and a priori probabilities in the route 2 MAP1.

The cESI estimator is then ${\overline{\Sigma }}_{\widehat{\iota }\widehat{\iota }}\left(f\right)$ (Eq. 21) for any set of inverse solution instances $\left({\widehat{\iota }}_m\left(f\right);\forall m\right)$, or more compactly from an estimate of the data cross-spectrum ${\overline{\Sigma }}_{vv}\left(f\right)$ (Eq. 22).

$\begin{array}{l}{\overline{\Sigma }}_{\widehat{\iota }\widehat{\iota }}\left(f\right)=\langle {\widehat{\iota }}_m\left(f\right){\widehat{\iota }}_m^{†}\left(f\right);M\rangle \langle {\widehat{\iota }}_m\left(f\right){\widehat{\iota }}_m^{†}\left(f\right);M\rangle =\\ \hfill \left(\frac{1}{M}\right){\displaystyle \sum _{m=1}^M}{\widehat{\iota }}_m\left(f\right){\widehat{\iota }}_m^{†}\left(f\right)& (21)\end{array}$

$\begin{array}{ll}{\overline{\Sigma }}_{\widehat{\iota }\widehat{\iota }}\left(f\right)={\widehat{\text{T}}}_{\iota v}\left(f\right){\overline{\Sigma }}_{vv}\left(f\right){\widehat{\text{T}}}_{v\iota }\left(f\right)& (22)\end{array}$

Important examples of cESI that follow route 1 are MNE, eLORETA and LCMV:

• The Minimum Norm Estimate (MNE) (Hämäläinen and Ilmoniemi, 1994). The basic smooth model of the source variables, that could either disregard the weight matrix A(f) = I or consider it ad-hoc A(f) = A^ac based on anatomical information.

• The Linearly Constrained Minimum Variance (LCMV) (Van Veen et al., 1997). The beamformer method that estimates a diagonal weight matrix A(f) = diag(a). This estimation produces an ideal filter (inverse operator) for each source variable, suppressing the interference from the other source variables and performing ideally under focalized distribution around one or a few sources.

• The Exact Low-Resolution Electromagnetic Tomographic Analysis (ELORETA) (Pascual-Marqui et al., 2006). A regression method that estimates A(f) so that the localization of the maximum for the estimated source variables corresponds exactly to the true maximum. This estimation performs ideally under a unimodal and smooth distribution of source variables.

The variants of these techniques are the one optimized for route 2, ESI in the frequency domain: the Spectral eLORETA (seLORETA) (Nolte et al., 2020) and the Spectral LCMV (sLCMV) (Larson-Prior et al., 2013). An additional solution used in this study as a reference is the spectral MNE (sMNE).

Novel cESI theory leading to the sSSBL approximation

Having described the theory of state-of-the-art cESI methods, we now focus on more sophisticated MAP theory and the sSSBL approximation that allows their practical implementation.

Bayesian MAP2 inverse operators

It is important to stress that posing ad-hoc priors and their hyperparameters is always necessary due to the uncertainties involved, by definition of the MAPs (Eqs. 14, 15) [54], [126]. In route 2, the prior was placed upon the Fourier transform ι(f) (Eq. 19). An alternative is to leverage the asymptotic distribution of the Fourier transform (Eq. 10) and place the a priori upon Σ_ιι(f).

This approach brings us to our main contribution, the theoretically promising cESI, which leads to an optimal inverse operator ${\widehat{\text{W}}}_{\iota v}$ based on the second-type MAP (MAP2) (Eq. 23).

$\begin{array}{ll}\begin{array}{c}\left(route3\right){\widehat{\Sigma }}_{\iota \iota }\left(f\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left({\overline{\Sigma }}_{vv}\left(f\right)\right)\\ {\widehat{\text{W}}}_{\iota v}=argma{x}_{{\text{W}}_{\iota v}}q\left({\text{W}}_{\iota v}\left({\overline{\Sigma }}_{vv}\left(f\right)\right)|{\overline{\Sigma }}_{vv}\left(f\right)\right)\\ q\left({\Sigma }_{\iota \iota }\left(f\right)|{\overline{\Sigma }}_{vv}\left(f\right)\right)\propto p\left({\Sigma }_{vv}\left(f\right)|{\Sigma }_{\iota \iota }\left(f\right)\right)p\left({\Sigma }_{\iota \iota }\left(f\right)\right)\end{array}& (23)\end{array}$

This MAP2 above posits the a priori probability p(Σ_ιι (f)) upon the source cross-spectrum Σ_ιι(f), considered a random hyper-parameter matrix. Incorporating an a priori to match the likelihood upon the sampled estimator ${\overline{\Sigma }}_{vv}\left(f\right)$ leads to a posteriori $q\left({\Sigma }_{\iota \iota }\left(f\right)|{\overline{\Sigma }}_{vv}\left(f\right)\right)$. Such a likelihood is the Complex Wishart probability density W^ℂ (Eq. 11) now defined for ${\overline{\Sigma }}_{vv}\left(f\right)$ (Eq. 24) with scale matrix Σ_vv(f) and M degrees of freedom. The specific scale matrix Σ_vv(f) in this likelihood posits the relation to the matrix equation of source cross-spectrum Σ_ιι(f) (Eq. 12). An a priori probability is then upon some cross-spectral Gibbs energy P(Σ_ιι (f)) defined by norms such as the vectorized (entry-wise) p-norm (Ding et al., 2006) and Schatten p-norm (Fan, 1951; Schatten, 2013) upon the matrix Σ_ιι(f).

$\begin{array}{ll}\begin{array}{c}p\left({\overline{\Sigma }}_{vv}\left(f\right)|{\Sigma }_{\iota \iota }\left(f\right)\right)=W^{ℂ}\left({\overline{\Sigma }}_{vv}\left(f\right)|{\text{L}}_{v\iota }{\Sigma }_{\iota \iota }\left(f\right){\text{L}}_{\iota v}+{\Sigma }_{\xi \xi }\left(f\right),M\right)\\ p\left({\Sigma }_{\iota \iota }\left(f\right)\right)\propto exp\left(H\left({\Sigma }_{\iota \iota }\left(f\right)\right)|\alpha \left(f\right)\right)\end{array}& (24)\end{array}$

Depending on the likelihood and a priori probability (Eq. 24), the optimal inverse operator ${\widehat{\text{W}}}_{\iota v}$ (Eq. 23) is often non-linear and numerically intractable. In this case, which we followed in this article, optimizing ${\widehat{\text{W}}}_{\iota v}$ requires Approximated Bayesian Computation (ABC) (Csilléry et al., 2010). The ABC we describe here employed ${\widehat{\text{W}}}_{\iota v}^{\left(k+1\right)}$, a non-convex but numerically tractable successive approximation to ${\widehat{\text{W}}}_{\iota v}$. In turn, the tractable ${\widehat{\text{W}}}_{\iota v}^{\left(k+1\right)}$ followed from $q^{\left(k\right)}\left({\Sigma }_{\iota \iota }\left(f\right)|{\overline{\Sigma }}_{vv}\left(f\right)\right)$, the non-convex relaxation of the a posteriori $q\left({\Sigma }_{\iota \iota }\left(f\right)|{\overline{\Sigma }}_{vv}\left(f\right)\right)$ in Expectation-Maximization (EM) iterations. Obtaining the a posteriori relaxation is via the Variational Bayes (VB) treatment (Dempster et al., 1977; Liu and Rubin, 1994; Daunizeau and Friston, 2007; Friston et al., 2007; Nummenmaa et al., 2007; Friston K. J. et al., 2008). Such a VB treatment could target the separable model for $q\left({\Sigma }_{\iota \iota }\left(f\right)|{\overline{\Sigma }}_{vv}\left(f\right)\right)$, or the separable Hierarchical Bayesian (HB) model for the a priori p(Σ_ιι (f)).

Bayesian (joint-MAP) inverse operators and cross-spectral norms

Considering the MAP2 (Eq. 23) and the latent cross-spectrum matrix Σ_ιι(f) constrained within the vector subspace generated by random instances (ι_m (f); ∀ m), i.e., as if it was defined by its usual estimator ${\Sigma }_{\iota \iota }\left(f\right)\cong {\overline{\Sigma }}_{\iota \iota }\left(f\right)$ (Eq. 7). Hence, the MAP2 is in probability equivalent to the joint-MAP (Hsiao et al., 1998, 2002; Yeredor, 2000; Davis et al., 2001; Auranen et al., 2005; Chen et al., 2011).

We introduce this joint-MAP (Eq. 25) substituting in the general MAP definition (Eqs. 14, 15) the data ${Y}$ and source variable ${X}$ by the instances ${Y}=\left(v_m\left(f\right);\forall m\right)$ and ${X}=\left({\iota }_m\left(f\right);\forall m\right)$.

$\begin{array}{ll}\begin{array}{c}\left(\widehat{{\iota }_m\left(f\right);\forall m}\right)\leftarrow {\widehat{\text{W}}}_{\iota v}\left(v_m\left(f\right);\forall m\right)\\ {\widehat{\text{W}}}_{\iota v}\left(v_m\left(f\right);\forall m\right)=argma{x}_{{\text{W}}_{\iota v}}p\left({\text{W}}_{\iota v}\left(v_m\left(f\right);\forall m\right)|v_m\left(f\right);\forall m\right)\\ p\left({\iota }_m\left(f\right);\forall m|v_m\left(f\right);\forall m\right)\propto p\left(v_m\left(f\right);\forall m|{\iota }_m\left(f\right);\forall m\right)p\left({\iota }_m\left(f\right);\forall m\right)\end{array}& (25)\end{array}$

In our case, we considered the above p(v_m(f);∀m|ι_m(f);∀m) the factorizable joint likelihood (Eq. 26) whose factors are present throughout the complex-valued Gaussian probability density (likelihood for route 2) (Eq. 18). In addition, one may assume Σ_ξξ = diag(β (f)) a univariate (in many cases homogenous) noise model as expressed in terms of the noise spectrum β(f).

$\begin{array}{ll}\begin{array}{c}p\left(v_m\left(f\right);\forall m|{\iota }_m\left(f\right);\forall m\right)={\displaystyle {\prod }_{m=1}^M}p\left(v_m\left(f\right)|{\iota }_m\left(f\right)\right)\\ p\left(v_m\left(f\right)|{\iota }_m\left(f\right)\right)=N^{ℂ}\left(v_m\left(f\right)|{\text{L}}_{v\iota }{\iota }_m\left(f\right),diag\left(\beta \left(f\right)\right)\right)\end{array}& (26)\end{array}$

We now introduce the joint a priori p(ι_m (f); ∀ m) (Eq. 27) upon the Gibbs energy H(ι_m (f); ∀ m) redefining the previous (Eq. 15) over instances (ι_m (f); ∀ m). Toward cESI, this function may be read as cross-spectral Gibbs energy $H\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)$.

$\begin{array}{ll}\begin{array}{c}p\left({\iota }_m\left(f\right);\forall m\right)\propto exp\left(H\left({\iota }_m\left(f\right);\forall m\right)|\alpha \left(f\right)\right)=\\ exp\left(H\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)|\alpha \left(f\right)\right)\end{array}& (27)\end{array}$

Our purpose now is to define the type of cross-spectral Gibbs energy that must deal with a severely ill-conditioned and high dimensional cESI setup (G ≫ S). We shall deal with these problems employing the vector or matrix norms, such as the vector structured p, q-norm (Kowalski and Torrésani, 2009a,b; Gramfort et al., 2012, 2013) and the Schatten matrix p-norm (Fan, 1951; Schatten, 2013).

To recap, our approach toward cESI is then with the joint-MAP inverse operator ${\widehat{\text{W}}}_{\iota v}$ (Eqs. 25–27) [instead of the MAP2 inverse operator (Eqs. 23, 24)]. Toward our target (cESI), we must leverage the class of joint a priori probabilities expressed by the cross-spectral Gibbs energy $H\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)$, not the more general case of Gibbs energy defined upon the instances H(ι_m (f); ∀ m).

We leverage H₁ (Eq. 28) the structured p = 1, q = 2-norm (square) that performs sparse regularization of the cross-spectrum topographic projection or spectrum $t{r}^{\frac{1}{2}}\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)$. H₁ is the well-known Least Absolute Shrinkage and Selection Operator (LASSO) (Tibshirani, 1996) for vectors ι(f) extended or structured over the second dimension index m for the set of instances (ι_m (f); ∀ m). Such a structured norm is a matrix quasinorm that does not fulfill the triangle equality over the cross-spectrum ${\overline{\Sigma }}_{\iota \iota }\left(f\right)$.

$\begin{array}{ll}{H}_1\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)={\sum }_{g=1}^G{\left({\sum }_{m=1}^M|{\iota }_m\left(g,f\right){|}^2\right)}^{\frac{1}{2}}=t{r}^{\frac{1}{2}}\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)& (28)\end{array}$

In addition, we introduce H₂ (Eq. 29) the structured p = 2, q = 2-norm (square), Shatten p = 1-norm or nuclear norm $tr\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)$, to compensate sparse bias of H₁ (Eq. 28) and regularize eigenvalues for the cross-spectrum ${\overline{\Sigma }}_{\iota \iota }\left(f\right)$. H₂ is the well-known Ridge operator (Hoerl and Kennard, 1970) for a vector ι(f) structured over the second dimension index m for the set of instances (ι_m (f); ∀ m).

$\begin{array}{ll}{H}_2\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)={\sum }_{m=1}^M{\sum }_{g=1}^G|{\iota }_m\left(g,f\right){|}^2=tr\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)& (29)\end{array}$

ESI practice, which is based on either the sparse or the smooth models, may be insufficient, though, in many scenarios where brain activity is patch-wise or not wholly sparse, Elastic Net, the linear combination of sparse/smooth models, may be ideal (Zou and Hastie, 2005; Vega-Hernández et al., 2008). Such a regularization style is a spatially structured sparsity due to the linear combination of the sparse LASSO and smooth Ridge operators upon the vector ι(f). Hence, for the vector instances (ι_m (f); ∀ m), the linear combination of the quasinorm H₁ (Eq. 28) and the nuclear norm H₂ (Eq. 29) leads to the following joint a priori probability p(ι_m (f); ∀ m) (Eq. 30).

$\begin{array}{ll}\begin{array}{c}p\left({\iota }_m\left(f\right);\forall m\right)\propto exp\left(H_1\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)|{\alpha }_1\left(f\right)\right)\\ exp\left(H_2\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)|{\alpha }_2\left(f\right)\right)\end{array}& (30)\end{array}$

Motivated by ESI practice, our definition of the Elastic Net (Eq. 30) (nuclear quasinorm) diverges from that of previous studies. From these works, the Elastic Net nuclear norm combines the nuclear norm $tr\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)$ with the square Schatten p = 2-norm (Frobenius norm) $tr\left({\overline{\Sigma }}_{\iota \iota }^2\left(f\right)\right)$ (Sun and Zhang, 2012; Chen et al., 2013; Kim et al., 2015; Zhang et al., 2017). This Elastic Net nuclear norm resolves convexity problems of the nuclear norm in the context of matrix completion (Candes and Recht, 2012; Hu et al., 2012) due to a non-convex problem with the sole nuclear norm. Although it was not our purpose to investigate convexity here, this property holds for our Elastic Net nuclear quasinorm, assuming cross-spectrum ${\overline{\Sigma }}_{\iota \iota }\left(f\right)$ upon vector basis (Zou and Hastie, 2005). In addition, we did not consider the vectorized (entry-wise) p-norm (Ding et al., 2006). The latter, which might be necessary to regularize off-diagonal entries in some cases (Paz-Linares et al., 2018), failed to ameliorate ill-condition or distortions.

The type of ABC introduced here is known as “gamma-MAP,” the standard VB treatment applied to a priori probabilities in the joint-MAP. In turn, the gamma-MAP leads to the quasilinear successive approximations ${\widehat{\text{T}}}_{\iota v}^{\left(k\right)}$. Solving the gamma-MAP under cross-spectral Gibbs energy $H\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)$ defined by the Elastic Net nuclear quasinorm is the Structured Sparse Bayesian learning algorithm (ssSBL) described in the next section.

Gamma-MAP and implementation of SSBL

Similarly to the MAP2 (Eq. 23), the joint-MAP inverse operator ${\widehat{\text{W}}}_{\iota v}$ (Eq. 25) could be non-linear by nature, depending on the joint a priori probability p(ι_m (f); ∀ m). We now introduce the gamma-MAP that achieves the quasilinear joint-MAP version leveraging the idea of mean-field approximation (Kadanoff, 2009). Such a mean-field approximation is the main idea behind the VB treatment applied to the latent variables of neuroimaging data (MacKay, 1998; Roweis and Ghahramani, 1999; Ghahramani and Beal, 2000, 2001; Eyink et al., 2004; Friston et al., 2007; Nummenmaa et al., 2007; Friston K. J. et al., 2008; Babacan et al., 2009).

By assuming the definition for this joint a priori p(ι_m (f); ∀ m) upon Gibbs energy H(ι_m (f); ∀ m) (Eq. 27), the mean-field approximation is then to search for the self-consistent energy form $\tilde{H}\left({\iota }_m\left(f\right);\forall m\right)$ (Eq. 31) as represented by separable (in most cases also indistinguishable) energy terms H_g(ι_m (g , f); ∀ m) for each source. A self-consistent energy form H_g(ι_m (g , f); ∀ m) also summarizes the interaction field g ↔ ∀g′, as a property of the original Gibbs energy H(ι_m (f); ∀ m).

$\begin{array}{ll}\begin{array}{c}H\left({\iota }_m\left(f\right);\forall m\right)\cong \tilde{H}\left({\iota }_m\left(f\right);\forall m\right)=\\ {\displaystyle {\sum }_g}{H}_g\left({\iota }_m\left(g,f\right);\forall m\right)\end{array}& (31)\end{array}$

Obtaining the self-consistent form Eq. (31) was always plausible under the Gibbs energy H(ι_m (f); ∀ m) that fulfills pair-wise separability in identical functions $H\left({\iota }_m\left(g,f\right),{\iota }_m\left(g^{{}^{\prime }},f\right);\forall m\right)$ (Eq. 32). Note that this property is known for ensuring the convergence of solutions for the self-consistency field equations in the literature of magnetism (Weiss, 1907, 2001; Le Boudec et al., 2007; Kadanoff, 2009; Zheng et al., 2015). The cross-spectral Gibbs energy $H\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)$, as defined by the norms considered here (Elastic nuclear quasinorm (Eq. 30) or others), fulfills such a property.

$\begin{array}{ll}\begin{array}{c}H\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)=\\ {\displaystyle {\sum }_{g=1}^G}{\displaystyle {\sum }_{g^{{}^{\prime }}=1}^G}H\left(\langle {\iota }_m\left(g,f\right){\iota }_m^{†}\left(g^{{}^{\prime }},f\right);M\rangle \right)\\ \end{array}& (32)\end{array}$

Then, if the joint a priori probability p(ι_m (f); ∀ m) (Eq. 27) is written in terms of the separable Gibbs energy $H\left({\overline{\Sigma }}_{\iota \iota }\left(f\right)\right)$ (Eq. 32) it turns out the Markov random field property (Eq. 33) that is summarized by $p\left({\iota }_m\left(g,f\right);\forall m|{\iota }_m\left(g^{{}^{\prime }},f\right);\forall m\right)$ the pair-wise and identical probability factors (Kindermann et al., 1980; Lafferty et al., 2001). Note that such probability factors may approximate but do not strictly represent conditional probabilities.

$\begin{array}{ll}\begin{array}{c}p\left({\iota }_m\left(f\right);\forall m\right)=\\ {\displaystyle {\prod }_{g=1}^G}{\displaystyle {\prod }_{g^{{}^{\prime }}=1}^G}p\left({\iota }_m\left(g,f\right);\forall m|{\iota }_m\left(g^{{}^{\prime }},f\right);\forall m\right)\\ p\left({\iota }_m\left(g,f\right);\forall m|{\iota }_m\left(g^{{}^{\prime }},f\right);\forall m\right)\propto \\ exp\left(P\left(\langle {\iota }_m\left(g,f\right){\iota }_m^{†}\left(g^{{}^{\prime }},f\right);M\rangle \right)|\alpha \left(f\right)\right)\\ \end{array}& (33)\end{array}$

Hence, we target VB—the factorizable (separable) approximation $h^q\left({\iota }_m\left(f\right);\forall m\right)$ (Eq. 34) of the joint a priori p(ι_m (f); ∀ m) (Eq. 33). Essential to obtain $h_g^q\left({\iota }_m\left(g,f\right);\forall m\right)$ is the Hierarchical Bayes (HB) or probability mixture bellow, identical a priori h and a posteriori q probabilities that are upon some type of “variational” hyper-parameters γ(g , f). Such an a posteriori q is also known as the proxy for belief propagation or message passing, denominated iterated conditional mode in the general literature of Markov random fields (Pearl, 1988, 2022; Weiss, 2001; Yedidia et al., 2003; Zheng et al., 2015).

$\begin{array}{ll}\begin{array}{c}p\left({\iota }_m\left(f\right);\forall m\right)\cong h^q\left({\iota }_m\left(f\right);\forall m\right)\propto {\displaystyle {\prod }_{g=1}^G}{h}_g^q\left({\iota }_m\left(g,f\right);\forall m\right)\\ h_g^q\left({\iota }_m\left(g,f\right);\forall m\right)=\\ \int h\left({\iota }_m\left(g,f\right);\forall m|\gamma \left(g,f\right)\right)q\left(\gamma \left(g,f\right)|{\iota }_m\left(f\right);\forall m\right)d\gamma \left(g,f\right)\\ \end{array}& (34)\end{array}$

Where the hyper-parameters γ(g , f) condense field interactions H(〈ι (g , f) ι^† (g′ , f); M〉) (Eq. 33) by the definition of an a posteriori probability q upon (ι_m (f); ∀ m).

After specifying an a priori h, the VB treatment is applied to search for estimators of the a posteriori $\widehat{q}$ (Eq. 35). That is to minimize the Kullback-Leibler divergence $D_{KL}\left(p,h^q\right)$ of the approximation h^qregarding to the joint a priori p (Eq. 34).

$\begin{array}{ll}\begin{array}{c}\widehat{q}=argmi{n}_q{D}_{KL}\left(p\left|\right|h^q\right)\\ D_{KL}\left(p\left|\right|h^q\right)=\\ \int p\left({\iota }_m\left(f\right);\forall m\right)log\left(p\left({\iota }_m\left(f\right);\forall m\right)/h^q\left({\iota }_m\left(f\right);\forall m\right)\right)\\ d\left({\iota }_m\left(f\right);\forall m\right)\end{array}& (35)\end{array}$

Minimizing the $D_{KL}\left(p,h^q\right)$, although theoretically achievable, turned out to be a difficult variational calculus problem for non-parametric HB (mixture) models (Blei and Jordan, 2006; Bryant and Sudderth, 2012; Gershman and Blei, 2012; Gershman et al., 2012; Nguyen and Bonilla, 2013; Duvenaud et al., 2016). Thus, the D_KL approach is common with parametric solutions constraining models (a priori h and a posteriori q) to the exponential family of probabilities (Andersen, 1970; Casella and Berger, 2021).

A simplifying assumption that could bypass the Kullback-Leibler divergence $D_{KL}\left(p,h^q\right)$ (Eq. 35) is to regard the a priori h and a posteriori q (34) within a family of conjugate HB models of the likelihood. Furthermore, we considered a priori h factorizable over the set of instances since the pair-wise interactions are additive by definition of 〈ι (g , f) ι^† (g′ , f); M〉 the cross-spectrum entries in Eq. (7).

Henceforth, a complex-valued Gaussian probability N^ℂ defines the a priori h(ι_m(f)|γ(f)), and a probability in the Gamma family Γ defines the a posteriori q(γ(f)|ι_m(f) ; ∀m) which we incorporate into the following HB model (Eq. 36).

$\begin{array}{ll}\begin{array}{c}{h}^q\left({\iota }_m\left(f\right);\forall m\right)\propto \int \left({\displaystyle {\prod }_{m=1}^M}h\left({\iota }_m\left(f\right)|\gamma \left(f\right)\right)\right)q\left(\gamma \left(f\right)|{\iota }_m\left(f\right);\forall m\right)\\ d\gamma \left(f\right)\\ h\left({\iota }_m\left(f\right)|\gamma \left(f\right)\right)=N^{ℂ}\left({\iota }_m\left(f\right)|0,diag\left(\gamma \left(f\right)\right)\right)\\ q\left(\gamma \left(f\right)|{\iota }_m\left(f\right);\forall m\right)={\displaystyle {\prod }_{g=1}^G}\Gamma \left(\gamma \left(g,f\right)|\delta \left(g,{\iota }_m\left(f\right);\forall m\right)\right)\end{array}& (36)\end{array}$

The previous HB model is the Generalized Gaussian Scale Mixture Model (GGSMM) prescribed by the Andrews and Mallows lemma (Andrews, 1974), where the particular and “Gamma” probability Γ depends on joint a priori to be approximated via $D_{KL}\left(p,h^q\right)$ (Eq. 35) and falls within the class of scale mixture Gamma models (McLachlan and Basford, 1988; Lindsay, 1995; Böhning and Seidel, 2003; Hancock and Samuelsen, 2007).

In such a mixture model, the hyper-parameter vector γ(f) may then be interpreted as the “variational spectrum,” which specifies a Complex Gaussian a priori probability N^ℂ (Eq. 36). The a posteriori for the variational spectrum entries γ(g , f) is a form of probability that belongs to the Gamma Γ family, with a parameterization δ(g , ι_m (f); ∀ m) that condenses field interactions H(〈ι (g , f) ι^† (g′ , f); M〉) (Eq. 32).

The joint-MAP (Eq. 25) can be approximated sequentially within iterations for the set $\left({\iota }_m^{\left(k\right)}\left(f\right);\forall m\right)$, or directly the cross-spectrum ${\widehat{\Sigma }}_{\iota \iota }^{\left(k\right)}\left(f\right)$, due to the quasilinear inverse operator ${\text{T}}_{\iota v}^{\left(k\right)}$ (37). Quasilinear ${\text{T}}_{\iota v}^{\left(k\right)}$ was then equivalent to iterated MAP1s for the Fourier transform (Eq. 20), positing complex-valued Gaussian likelihood p(v_m(f)|ι_m(f)) (Eq. 26) and a priori h(ι_m(f)|γ(f)) (36).

$\begin{array}{ll}\begin{array}{c}\widehat{\left({\iota }_m^{\left(k\right)}\left(f\right);\forall m\right)}\leftarrow {\text{T}}_{\iota v}^{\left(k\right)}\left(f\right)\left(v_m\left(f\right);\forall m\right)\\ {\widehat{\overline{\Sigma }}}_{\iota \iota }^{\left(k\right)}\left(f\right)={\text{T}}_{\iota v}^{\left(k\right)}\left(f\right){\overline{\Sigma }}_{vv}\left(f\right){\text{T}}_{v\iota }^{\left(k\right)}\left(f\right)\\ \end{array}& (37)\end{array}$

These iterated MAP1s are for an equivalent a posteriori $q^{\left(k\right)}\left({\iota }_m\left(f\right)|v_m\left(f\right)\right)$ (Eq. 38), specified by the mean ${\text{T}}_{\iota v}^{\left(k\right)}\left(f\right)v\left(f\right)$ and the covariance ${\Pi }_{\iota \iota }^{\left(k\right)}\left(f\right)$, which are functions of the variational spectrum diag(γ^(k) (f)) and the noise spectrum diag(β^(k) (f)). These spectrums specified the type of univariate approximations for the source cross-spectrum Σ_ιι , and the noise cross-spectrum Σ_ξξ in the quasilinear MAP1 of the Fourier transform (Eq. 20). Henceforth, we place the emphasis in the variational spectrum iterations γ^(k)(f) and defer this noise spectrum β^(k)(f) which is not essential for our main exposition.

$\begin{array}{ll}\begin{array}{c}{q}^{\left(k\right)}\left({\iota }_m\left(f\right)|v_m\left(f\right)\right)=N^{ℂ}\left({\iota }_m\left(f\right)|{\text{T}}_{\iota v}^{\left(k\right)}\left(f\right)v_m\left(f\right),{\Pi }_{\iota \iota }^{\left(k\right)}\left(f\right)\right)\\ {\text{T}}_{\iota v}^{\left(k\right)}\left(f\right)={\widehat{\Pi }}_{\iota \iota }^{\left(k\right)}\left(f\right){\text{L}}_{\iota v}dia{g}^{-1}\left({\beta }^{\left(k\right)}\left(f\right)\right)\\ {\Pi }_{\iota \iota }^{\left(k\right)}\left(f\right)={\left({\text{L}}_{\iota v}dia{g}^{-1}\left({\beta }^{\left(k\right)}\left(f\right)\right){\text{L}}_{v\iota }+dia{g}^{-1}\left({\gamma }^{\left(k\right)}\left(f\right)\right)\right)}^{-1}\\ \end{array}& (38)\end{array}$

Targeting the variational vector γ(f) was under the marginal or hyper-parametrized likelihood for data p(v_m(f)|γ(f)) (Eq. 39). This likelihood was due to integration (expectation) of p(v_m(f)|ι_m(f)) under the a priori probability for parameters h(ι_m(f)|γ(f)), which is upon the variational hyper-parameters (spectrum) γ(f) (Eq. 30) (MacKay, 1999).

$\begin{array}{ll}\begin{array}{c}p\left(v_m\left(f\right)|\gamma \left(f\right)\right)=\int p\left(v_m\left(f\right)|{\iota }_m\left(f\right)\right)h\left({\iota }_m\left(f\right)|\gamma \left(f\right)\right)d{\iota }_m\left(f\right)\\ p\left(v_m\left(f\right)|{\iota }_m\left(f\right)\right)=N^{ℂ}\left(v_m\left(f\right)|{\text{L}}_{v\iota }{\iota }_m\left(f\right),diag\left(\beta \left(f\right)\right)\right)\\ h\left({\iota }_m\left(f\right)|\gamma \left(f\right)\right)=N^{ℂ}\left({\iota }_m\left(f\right)|0,diag\left(\gamma \left(f\right)\right)\right)\\ \end{array}& (39)\end{array}$

However, the actual marginal likelihood p(v_m(f)|γ(f)) (Eq. 39) was intractable, with approximations via the expectation of the log joint probability p(v_m(f), ι_m(f)|γ(f)) (Dempster et al., 1977; Liu and Rubin, 1994) under the iterated a posteriori $q^{\left(k\right)}\left({\iota }_m\left(f\right)|v_m\left(f\right)\right)$ (Eq. 40). Then an approximated marginal likelihood $L\left(v_m\left(f\right)|{\gamma }^{\left(k\right)}\left(f\right),\gamma \left(f\right)\right)$ depended on the variational hyper-parameters (spectrum) γ^(k)(f)in the current iteration.

$\begin{array}{l}\begin{array}{c}log\left(p\left(v_m\left(f\right)|\gamma \left(f\right)\right)\right)\cong log\left(L\left(v_m\left(f\right)|{\gamma }^{\left(k\right)}\left(f\right),\gamma \left(f\right)\right)\right)=\\ \int q^{\left(k\right)}\left({\iota }_m\left(f\right)|v_m\left(f\right)\right)log\left(p\left(v_m\left(f\right)|{\text{L}}_{v\iota }{\iota }_m\left(f\right)\right)h\left({\iota }_m\left(f\right)|\gamma \left(f\right)\right)\right)\end{array}\\ \hfill d{\iota }_m\left(f\right)\hfill & (40)\end{array}$

An inverse operator W_ιv of the variational spectrum γ(f) is theoretically equivalent to the previous joint-MAP (Eq. 25), which is inverse operator of the samples (ι_m (f); ∀ m). This is known from the literature as “gamma-MAP” (Hsiao et al., 1998, 2002; Wipf and Nagarajan, 2009). Within the loop (Eq. 41) such an inverse operator ${\text{W}}_{\iota v}^{\left(k\right)}$ was the tractable successive approximations, with the iterated a posteriori $q\left(\gamma \left(f\right)|{\iota }_m^{\left(k\right)}\left(f\right);\forall m\right)$ and likelihood $L\left(v_m\left(f\right)|{\gamma }^{\left(k\right)}\left(f\right),\gamma \left(f\right)\right)$.

$\begin{array}{ll}\begin{array}{c}{\gamma }^{\left(k+1\right)}\left(f\right)\leftarrow {\widehat{\text{W}}}_{\iota v}^{\left(k\right)}\left(v_m\left(f\right);\forall m\right)\\ {\widehat{\text{W}}}_{\iota v}\left(v_m\left(f\right);\forall m\right)=\\ argma{x}_{{\text{W}}_{\iota v}}{p}^{\left(k\right)}\left({\text{W}}_{\iota v}\left(v_m\left(f\right);\forall m\right)|v_m\left(f\right);\forall m\right)\\ p^{\left(k\right)}\left(\gamma \left(f\right)|v_m\left(f\right);\forall m\right)\propto \left({\displaystyle {\prod }_{m=1}^M}L\left(v_m\left(f\right)|{\gamma }^{\left(k\right)}\left(f\right),\gamma \left(f\right)\right)\right)\\ q\left(\gamma \left(f\right)|{\iota }_m^{\left(k\right)}\left(f\right);\forall m\right)\end{array}& (41)\end{array}$

Here we implemented the Bayesian learning algorithm that effectuates the gamma-MAP γ^{(k + 1)}(f) (41) and computes the quasilinear ${\widehat{\text{T}}}_{\iota v}^{\left(k\right)}$ (Eq. 38), for a joint a priori based on the Elastic Net nuclear quasinorm (Eq. 30). For this joint a priori, the solution to the HB (mixture) model (Eq. 36) is exact, following an extension of the Andrews and Mallows lemma to the structured sparsity norms which are here upon the spectrum ${\overline{\sigma }}_{\iota \iota }^2\left(f\right)$.

Henceforth, we refer to this algorithm as Spectral Structured Sparse Bayesian Learning (ssSBL). The full derivation of the gamma-MAP and the ssSBL algorithm are developed in Supplementary material (Section Standard cESI theory).

Comparison of the distortions produced by seLORETA, sLCMV, and sSSBL

Data used for the validation

We used two different sets of EEG data to calculate estimated CST, each with their corresponding gold standard CST:

(1) Realistically simulated low-density EEG obtained considering as sources a CST obtained from a MEG recording with a very high sensor density (simulated-EEGvsMEG). The gold standard here is the sMNE CST of the MEG recording due to the well-known advantages of MEG over EEG and the very high sensor density.

(2) Real low-density macaque EEG (EEGvsECoG). The gold standard here is the sMNE CST of the simultaneously recorded EcoG.

Simulated-EEG vs. MEG

Figure 3 was based on a high-quality resting state MEG recording for subject 175,237 from the HCP database. The MEG signal selected for this purpose (Figure 3a2) was the 246-channel preprocessed data file. The electrical and magnetic lead fields were calculated with the subject's cortical and head geometry. With the magnetic lead field, Spectral MNE (sMNE) was used to calculate the MEG sources, which were taken as sources for the EEG. These sources were passed through the electrical Lead Field (Figure 3a3) to simulate a low-density 19-channel EEG recording (Figure 3a1). The simulation design is standard, essentially the same as those based on more straightforward configurations, where dipoles or patches are taken as the “ground truth” (Haufe and Ewald, 2019). Here, we used a much more realistic set of sources, determined from the MEG. Code availability: https://github.com/CCC-members/MEGvsSimulated-EEG.

FIGURE 3

Figure 3. Simulated-EEG vs. MEG experiment that illustrates our validation strategy. We show distortion measures of the estimated Cortical Spectral Topographies (CST) resting-state MEG alpha activity. The same strategy is followed for all spectral bands. It is also applied to the study of distortion based on the EEG vs. ECoG experiment. We start from the cross-spectrum data shown in hot colormaps of 2D space topographies (a1) ${\sum }_{vv}^{EEG}$ for EEG sensors (green dots) and (a2) ${\sum }_{vv}^{MEG}({\sum }_{vv}^{ECoG}$) for MEG (ECoG) sensors (blue dots). Using the corresponding cross-spectrum data and Lead Fields (a3) ${\text{L}}_{v\iota }^{EEG}$ and ${\text{L}}_{v\iota }^{MEG}$ (${\text{L}}_{v\iota }^{ECoG}$), for EEG sensors (green dots) and MEG (ECoG) sensors (blue dots), upon human (macaque) cortex, head layers geometries, we obtain the inverse operators for EEG ${\stackrel{\wedge }{T}}_M^{EEG}$, based on the tested method “M,” and MEG (ECoG) ${\stackrel{\wedge }{T}}_{sMNE}^{MEG}\left({\stackrel{\wedge }{T}}_{sMNE}^{ECoG}\right)$, based on the reference method sMNE. Employing these inverse operators, we obtain estimators for (b1) ${\sigma }_M^{EEG}$ the EEG tested CST given by method “M” and (b2) σ₀ the MEG (ECoG) sMNE reference CST. Incongruence between (b1) the EEG tested CST and (b2) the MEG (ECoG) reference CST is measured through (b3) e_M the Earth Movers' Distance (EMD) and c_M the correlation distance (1-CORR). Using (a1) the EEG cross-spectrum data and (b3) Lead Field for EEG, we obtain the resolution operator R_M. Leakage in (b1) the EEG tested CST, which is based in (c2) the 1st quartile point of (b2) the MEG (ECoG) reference CST, is measured through (c1) r_M, the Generalized Point Spread Function (GPSF) and b_M the Blurring for the GPSF.

EEG vs. ECoG

Figure 4 was based on high-density macaque ECoG recordings acquired in 128 sensors, concurrently with a low-density EEG acquired at 19 sensors simultaneously with the ECoG in the resting state (Nagasaka et al., 2011). This macaque preparation allowed the realistic measurement of distortions in resting-state connectivity estimated from low-density EEG using different aESI solutions (Wang et al., 2019). A sensor array placed surgically on the left macaque's cortical surface (Figure 4a1) allowed dense recordings of ECoG (Figure 4a2). The lead fields for EEG and ECoG (Figure 4a3) were obtained using the macaque's cortical and head geometry. The ECoG lead field was used to generate a reference or ground truth CST using sMNE. It was unnecessary to simulate data since the concurrent low-density EEG was available. Code availability: https://github.com/CCC-members/ECoGvsEEG.

FIGURE 4

Figure 4. Confirmation of Cortical Spectral Topographies (CST) based in EEG recorded simultaneously with (a1) ECoG implanted in the macaque. An X-ray image shows (a2) the high-density ECoG array placed onto the macaque cortical surface. ECoG recordings and (a3) their Lead Fields provide a more fine-grained reference for confirming CST estimators and measures of distortions for the EEG. The validation here includes elements analogous to the MEGvsSimulated-EEG experiment of Figure 3 (a3, b1–b3, c1, c2). Incongruence between (b1) σ_M the EEG tested CST and (b2) σ_o the ECoG reference CST is measured through (b3) $e_M^{ECoG\to EEG}$ the Earth Movers' Distance (EMD) and $C_M^{ECoG\to EEG}$ the correlation distance (1-CORR). Leakage in (b1) σ_M the EEG tested CST, which is based in (c2) on the first quartile point of (b2) σ_o the ECoG reference CST is measured through (c1) r_M, the Generalized Point Spread Function (GPSF) and b_M the Blurring for the GPSF. Elements (a1, a2) of this figure, are freely available in http://www.www.neurotycho.org/.

Validation rationale

We evaluated the low-density EEG recordings using ssSBL, sELORETA, and sLCMV to produce the CSTs σ_sSSBL, σ_sELOR, and σ_sLCMV, respectively. We denote these CST generically as σ_M, with Mspecifying the method. These calculations were carried out for the previously mentioned experiments: simulated EEG vs. MEG and EEG vs. ECOG. These CSTs were compared with reference CSTs (σ₀) considered the “gold standards.” This reference was varied according to the type of measurement and will be detailed in the context below. Finally, each method's distortion of its reference was assessed using the measures described next.

Measures of distortion

Measures of distortion to compare σ_M with σ₀ fall into two groups, leakage and incongruence indices that we illustrate in the Simulated-EEG vs. MEG experiment (Figure 3). The measures of the macaque EEG vs. ECoG experiment in figure elements (Figures 4b1–b3, c1, c2) are analogous to those of the Simulated-EEG vs. MEG experiment.

Leakage (Palva et al., 2018; Van de Steen et al., 2019), or spread, is quantified using the Generalized Point Spread Function and a blurring measure. These are based on the concept of the resolution operator R_M = T_{ιv, M}L_vι. Here T_{ιv, M} is the inverse operator for the method M, and L_vι the lead field. Consider a point source u_o indexed by g₀, any column R_M(:, g₀) is then the voltages v₀ produced by this source.

• The Generalized Point Spread Function (Grova et al., 2006b; Haufe et al., 2008) (GPSF) r_M in Figure 3c1, depicted with a hot colormap, represents the leakage (spread) of the low-density EEG CST estimators for a given set of cortical points. These cortical points, shown as blue dots (Figure 3c2), are selected as the most active 25% of the reference σ₀ and are the set G₀. The GPSF, for any point g ∈ G in the set of cortical points, is then:

$\begin{array}{ll}{\text{r}}_M\left(g\right)=\sqrt{\frac{1}{\left|{𝔾}_0\right|}{\sum }_{g_0\in {𝔾}_0}|{\text{R}}_M\left(g,g_0\right){|}^2}& (42)\end{array}$

Where |𝔾₀| denotes the number of elements in that set.

• The blurring measure for images (BLUR) b_M is defined as the Spatial Dispersion (SD) of the GPSF (Grova et al., 2006a; Haufe et al., 2008). It is worth clarifying that for a perfect cESI solution, with zero b_M = 0 in Figure 3c1, there would be an exact coincidence between the GPSF r_M non-zero values in the colormap and the blue dots in Figure 3c2. Formally, b_M is defined as:

$\begin{array}{ll}{b}_M=\sqrt{\frac{1}{\left|{𝔾}_0\right|}{\sum }_{g_0\in {𝔾}_0}{\vartheta }_M^2\left(g_0\right)},& (43)\end{array}$

Where ${\vartheta }_M^2\left(g_0\right)$ is the spatial dispersion around the reference point g₀ and

$\begin{array}{ll}{\vartheta }_M^2\left(g_0\right)=\frac{1}{\left|𝔾\right|}{\sum }_{g\in 𝔾}\left|{\text{R}}_M\left(g,g_0\right)\right|d_{g{g}_0}^2,& (44)\end{array}$

Where $d_{g{g}_0}^2$ is the square of the geodesic distance between those points.

Incongruence (Wang et al., 2019) quantifies the global level of distortion (leakage and localization error):

(1) The Earth Movers' Distance (EMD) e_M in Figure 3b3, which measures the effort to deform the CST spatial density determined from EEG (Figure 3b1) into the reference (Figure 3b2) (Grova et al., 2006b; Haufe et al., 2008; Paz-Linares et al., 2017).

(2) The correlation distance (1-CORR) c_M which measures deformations from the expected collinearity between pairs of spatially distributed CSTs.

These incongruence measures (EMD) e_M correlation distance (1-CORR) c_M combine sensitivity to leakage and localization error.

Results

Simulated-EEG vs. MEG inverse solutions

Both sELORETA and sLCMV CST estimators were seriously affected by leakage (Figure 5), judging by the mismatch between their estimated GPSF r_sELOR and r_sLCMV compared to the reference points (blue dots). They were centered at incorrect sites (opposite the blue dots) and with a much larger spread. For ssSBL, the set of local maxima for r_ssSBL was correctly centered around reference points and did not extend significantly beyond these. In other words, as can be appreciated qualitatively, r_ssSBL minimized leakage compared to sELORETA and LCMV.

FIGURE 5

Figure 5. In hot colormaps, measurements of “Leakage” in the EEG-based Cortical Spectral Topographies (CST) were obtained with all tested methods “M” (ssSBL, seLORETA, and sLCMV). Leakage is measured by employing r_M, the Generalized Point Spread Function (GPSF) regarding to the MEG sMNE reference CST, shown in the cortical views Left (L), Dorsal (D), Posterior (P), and Ventral (V) and for five characteristic spectral bands (delta, theta, alpha, beta, and gamma). Leakage distortions are proportional to mismatch and spread in the GPSF values (hot colormaps) regarding reference points (blue dots). Calculations of the GPSF follow the procedure described in Figure 3 for the MEGvsSimulated-EEG experiment.

Initially, these results differed from those reported by other authors for r_sELOR and r_sLCMV (Haufe and Ewald, 2019). The explanation may be due to their using idealized EEG simulations (point generators or discrete patches) that generate fewer distortions, thus reducing the apparent differences in the performance of different aESI methods. We verified this in recent aESI studies (Paz-Linares et al., 2017), concluding that simple SSBL (not the cESI ssSBL implemented in this paper) outperforms eLORETA and other methods only by a narrow margin.

An aESI solution computed with LCMV was usually sparser than an ELORETA solution, which was expected due to the thresholding strategy implemented in the original LCMV (not sLCMV) (Van Veen et al., 1997). However, this did not lead to any improvement in terms of the leakage observed in the GPSFs r_sELOR and r_sLCMV. The “zero localization error” of eLORETA has been claimed to be the unique advantage in its favor (Pascual-Marqui et al., 2006). However, this property has been theoretically demonstrated only for the peak (maximum) activity and not for the numerous local maxima of secondary activations common in real-life scenarios. Therefore, unsurprisingly, eLORETA can produce much better results in idealized simulations that use activity modeled as a focalized concentrated unimodal distribution (e.g., Gaussian) or only a few simple point sources (Kobayashi et al., 2003; Grova et al., 2008; Schoffelen and Gross, 2009; Haufe et al., 2013; Burle et al., 2015; Colclough et al., 2015; Bradley et al., 2016; Mahjoory et al., 2017; Stokes and Purdon, 2017; Palva et al., 2018; Haufe and Ewald, 2019; Marinazzo et al., 2019).

A striking observation was that despite σ₀ being highly frequency dependent, the GPSF pattern was r_seLORETA and r_sLCMV was relatively invariant for all frequencies across the spectrum. This invariance suggests that the pattern was mainly due to Lead Field properties rather than physiologically (Lopes da Silva et al., 1974; Niedermeyer and da Silva, 2005; Lopes da Silva, 2013). In contrast, the GPSF r_ssSBL pattern was closely tied to the physiological fluctuations across frequencies: from activity in the slow delta band to the faster alpha band. These fluctuations were inherited from the reference MEG data.

The quantitative analysis of the leakage measure b_M (BLUR) (Figure 6, left-column) confirms our qualitative impressions based on the GPSF r_M. The radar graphs showed a decrease in leakage of ssSBL compared to sELORETA and sLCMV (b_sELOR > b_ssSBL, b_sLCMV > b_ssSBL). This improvement was valid for all spectral bands (delta, theta, alpha, gamma 1, and gamma 2).

FIGURE 6

Figure 6. In radar graphs, measurements of “Leakage” and “Incongruence” in the Cortical Spectral Topographies (CST) were obtained with all tested methods “M” (ssSBL, seLORETA, and sLCMV) from the EEG. Leakage, regarding “Y,” the MEG or ECoG sMNE reference CST, is measured employing B_M, the Blurring (BLUR), and Incongruence, employing e_M the Earth Movers' Distance (EMD), and C_M the Correlation Distance (1-CORR), shown for five characteristic spectral bands (delta, theta, alpha, beta, and gamma). Distortions based in the BLUR, EMD, and 1-CORR are proportional to their values in the radar graph. Calculations of the BLUR, EMD, and 1-CORR follow the procedure described in Figures 3, 4, for the MEGvsSimulated-EEG and ECoGvsEEG experiments.

EMD values e_M are intuitively the amount of work to deform the EEG-based CST estimator to the reference CST estimator (σ_M → σ₀). They were one to two orders larger for sELORETA and sLCMV than SSBL (e_sELOR >> e_ssSBL, and e_sLCMV >> e_ssSBL).

A linear model adjusted to every source of the EEG and MEG CSTs estimated data pairs shows a clear linear tendency (Figure 7), with correlation distances larger for sELORETA and sLCMV. Some correlations are even negative correlations. By contrast, c_sELOR was in the range of around 0.7, for all frequency bands. This behavior was congruent with that observed in the colormaps for GPSF (r_sELOR and r_sLCMV) (Figure 5), where the maximum values appear in opposite areas to the blue dots (reference estimation).

FIGURE 7

Figure 7. Linear model and correlations for the Cortical Spectral Topographies (CST) obtained from the MEGvsSimulated-EEG experiment in the human. These were adjusted from the CST data pairs (σ_M, σ₀) for all tested methods “M” (ssSBL, seLORETA, and sLCMV), and in five characteristic spectral bands (delta, theta, alpha, beta, and gamma).

These results suggest that commonly used ESI validation procedures, limited to idealized simulations of local neural currents, may not accurately assess the actual distortions (Haufe et al., 2013; Haufe and Ewald, 2019). Our simulation of EEG that incorporates realistic local neural currents, estimated from high-density MEG, shows that the leakage and localization error in ESI applied to actual data might be much more severe than expected (Palva et al., 2018; He et al., 2019; Van de Steen et al., 2019). Therefore, future efforts should consider validation benchmarks based on realistic simulations like those described here. When using this benchmark, the ssSBL approach appears to considerably control the effect of distortions with respect to other techniques.

EEG vs. ECoG inverse solutions

The GPSF colormaps (Figure 8) for each method showed a consistent behavior as those of the Simulated-EEG vs. MEG experiment (Figure 5). As evident in the GPSF maps r_ssSBL and measured in b_ssSBL (Figure 6, right-column), the performance of ssSBL in reducing leakage was superior compared to sELORETA (b_sELOR > b_ssSBL) and sLCMV (b_sLCMV > b_ssSBL).

FIGURE 8

Figure 8. In hot colormaps, measurements of “Leakage” in the EEG-based Cortical Spectral Topographies (CST) were obtained with all tested methods “M” (ssSBL, seLORETA, and sLCMV). Leakage is measured by employing r_M, the Generalized Point Spread Function (GPSF) regarding to the ECoG sMNE reference CST, shown in the cortical views Left (L), Dorsal (D), Posterior (P), and Ventral (V) and for five characteristic spectral bands (delta, theta, alpha, beta, and gamma). Leakage distortions are proportional to mismatch and spread in the GPSF values (hot colormaps) regarding reference points (blue dots). Calculations of the GPSF follow the procedure described in Figure 4 for the ECoGvsEEG experiments.

The measurements of incongruence by the EMD e_M (Figure 6, right-column) were consistent with those of the previous experiment (left-column), confirming the improvement of ssSBL regarding sELORETA (e_sELOR > e_ssSBL by a considerably narrow margin) and LCMV (e_sLCMV >> e_ssSBL by a wider margin of three orders of magnitude) for all spectral bands.

The correlation distance c_M for all methods showed that linear regression described the relation of EEG CST to ECoG CST well. The correlation was positive for ssSBL but negative for sELORETA and sLCMV as seen in Figure 6 (right-column). As a consequence (c_sLCMV > c_LORETA>c_ssSBL) for all spectral bands. This linear tendency in Figure 9, similarly to Figure 7 confirms the feasibility of cESI, even with low-density EEG, given its close relationship to a more direct observation modality like ECoG. The results of sELORETA and sLCMV did not reveal a clear linear positive tendency. These results confirm and extend the results of previous studies (Marinazzo et al., 2019; Papadopoulou et al., 2019; Wang et al., 2019), suggesting that some types of ESI should be interpreted with extreme care.

FIGURE 9

Figure 9. Linear model and correlations for the Cortical Spectral Topographies (CST) obtained from the ECoGvsEEG experiment in the macaque. These were adjusted from the CST data pairs (σ_M, σ_M) for all tested methods “M” (ssSBL, seLORETA, and sLCMV), and in five characteristic spectral bands (delta, theta, alpha, beta, and gamma).

Discussion

We now summarize and evaluate our results from a theoretical point of view. We introduce a general Bayesian framework for cESI, the estimation of source cross-spectral matrices. This approach allowed us to address the high level of topographic distortions, which arise from the severely ill-conditioned nature of the underlying inverse problem (Kobayashi et al., 2003; Grova et al., 2008; Schoffelen and Gross, 2009; Haufe et al., 2013; Burle et al., 2015; Colclough et al., 2015; Bradley et al., 2016; Mahjoory et al., 2017; Stokes and Purdon, 2017; He et al., 2018, 2019; Palva et al., 2018; Haufe and Ewald, 2019; Marinazzo et al., 2019). We indicated that the problems originating from the ill-posedness were compounded by the habitual implementation of linear and non-sparse (smooth) types of inverse solutions (Mantini et al., 2007; Marzetti et al., 2008; Brookes et al., 2011b, 2012; Hipp et al., 2012; Lopes da Silva, 2013; Colclough et al., 2015; Marinazzo et al., 2019; Nolte et al., 2020).

Non-linear and sparse inverse solutions might ameliorate distortions (Tibshirani, 1996; Zou and Hastie, 2005; Yuan and Lin, 2006; Haufe et al., 2008; Vega-Hernández et al., 2008; Kowalski and Torrésani, 2009a,b; Li and Tian, 2011; Gramfort et al., 2012, 2013). Unfortunately, they may also introduce warping and other biases of the cross-spectral (cESI) estimator. ESI practice with these inverse solutions was most beneficial when addressing the deterministic time/frequency waveforms of event-related brain processes within the framework of spatially short-scale distributed event-related brain networks (Picton and Hillyard, 1974; Picton et al., 1974; Lopes da Silva et al., 1991; Clark et al., 1994; Makeig et al., 1999, 2004; Makeig, 2002; Eichele et al., 2005; Harrison et al., 2008; Vega-Hernández et al., 2008; Maurer and Dierks, 2012). This approach to ESI empiricism suggests that the best results were obtained with flexible smooth/sparse a priori models (Zou and Hastie, 2005; Vega-Hernández et al., 2008; Li and Lin, 2010). Quasilinear inverse solutions of these smooth/sparse models can provide good data-driven approximations.

We formalized the three possible routes toward cESI (our target) via MAP1 inverse solutions, for the vector processes or their Fourier transform and via MAP2 or joint-MAP inverse solutions for their cross-spectrum. The cross-spectral MAP2 or joint-MAP are plausible inverse solutions, which target the quantity of interest and theretofore posit the a priori upon the cross-spectrum Σ_ιι(f).

In contrast, most cESI have been previously limited to route 1 for the processes or route 2 for the Fourier transform (aESI) wherein the a priori is placed upon ι(t) or ι(f). Therefore, essential statistics like cross-spectrum were incorrectly addressed as the subsequent step to aESI. As we demonstrated in simulations and real data, the aESI procedure is not suitable, with cESI amplifying the distortions previously produced during aESI. Indeed, a way to achieve statistical guarantees in cESI is through route 3 associated with the source cross-spectrum for the data.

We have deferred for now to the complete implementation of the MAP2 inverse solution for the cross-spectrum. This MAP2, which targets a source matrix, was not straightforward, given a cESI setup that is severely ill-conditioned and of high dimensionality. Hence, we adopted the joint-MAP, as an approximation to the MAP2, which targets ${\overline{\Sigma }}_{\iota \iota }\left(f\right)$ of the sampled source cross-spectrum.

An essential concept introducing any MAP inverse solution was the Gibbs energy, a concept first formulated as the Gibbs free energy of any physical system (Landau and Lifshitz, 1980). From the point of view of the inverse problem theory, the Gibbs free energy describes the forward energy exchange of a system (the source variable) to the media (the data) (Ghahramani and Beal, 2001; Grave de Peralta Menendez et al., 2004; Friston et al., 2006; Mattout et al., 2006; Friston K. J. et al., 2008; Friston K. L. et al., 2008; Wipf and Nagarajan, 2009). Thus, this is an important generalization of the Tikhonov regularization and ESI for source cross-spectrum.

We employed this cross-spectral Gibbs energy to model the joint a priori probabilities and second-order multivariate properties of the Fourier transform. The Gibbs energy in cESI must be a function of the cross-spectrum entries. This assumption follows from the Gaussianity of Fourier transform (Brillinger, 1965; Brillinger and Rosenblatt, 1967), or statistical sufficiency of the cross-spectrum. This assumption is valid for ESI under a variety of experimental conditions, as it can be demonstrated with high-quality MEG and ECoG data (Nagasaka et al., 2011; Van Essen et al., 2013).

Quasilinear inverse solutions preserved the F-invariance for Gaussianity, avoiding warping of cross-spectral amplitude and phase information (Marzetti et al., 2008; He et al., 2019; Nolte et al., 2020). F-invariance is also valid for a MAP1 assuming Gaussianity of the Fourier transform, and streamlining the dimensionality reduction for MAP2, leading to our joint-MAP interpretation (Hsiao et al., 1998, 2002; Yeredor, 2000; Davis et al., 2001; Auranen et al., 2005; Chen et al., 2011).

We indicate the importance of the nuclear norm (trace) and nuclear quasinorm (square root trace) for matrices, from the context of matrix completion inverse problems (Fan, 1951; Ding et al., 2006; Sun and Zhang, 2012; Chen et al., 2013; Schatten, 2013; Kim et al., 2015; Zhang et al., 2017) which translated into a sparse/smooth spectrum model. Per the definition of cross-spectrum upon vector basis, this model indeed could be unified with the structured vector norms (Kowalski and Torrésani, 2009a,b; Gramfort et al., 2012, 2013). Furthermore, this has been very common in aESI with sparse (LASSO) (Tibshirani, 1996) and smooth (MNE) (Hoerl and Kennard, 1970) model. Our sparse-smooth model was, therefore, an Elastic Net nuclear quasinorm.

Our approach has an important connection to Bayesian learning (Tipping, 2001; Wipf et al., 2006; Park and Casella, 2008; Wipf and Nagarajan, 2009; Casella et al., 2010; Li and Lin, 2010; Li et al., 2011; Paz-Linares et al., 2017) and provided the link between the joint-MAP and quasilinear inverse solutions. We introduced a type of variational approximation to the joint a priori probability upon cross-spectral Gibbs energy. This variational approximation is similar to the previous Bayesian Learning methods with extended applicability to high dimensional inverse problems that are solved in iterated conditional mode (Pearl, 1988, 2022; Weiss, 2001; Yedidia et al., 2003; Zheng et al., 2015).

Here, we restricted ourselves to the Gaussian Mixture Model approach (Blei and Jordan, 2006; Bryant and Sudderth, 2012; Gershman and Blei, 2012; Gershman et al., 2012; Nguyen and Bonilla, 2013; Duvenaud et al., 2016) applied to the specific a priori per the definition of the Elastic Net nuclear quasinorm (Candes and Recht, 2012; Hu et al., 2012). Our Bayesian learning method (ssSBL) was a generalization of the SBL and sSBL approaches.

Studying the distortions in cESI required a validation based on a reference estimation (ground truth) closer to actual source distributions in the brain rather than simulations of brain electrical activity based on ideal source configurations (Kobayashi et al., 2003; Grova et al., 2008; Schoffelen and Gross, 2009; Haufe et al., 2013; Burle et al., 2015; Colclough et al., 2015; Bradley et al., 2016; Mahjoory et al., 2017; Stokes and Purdon, 2017; Palva et al., 2018; Haufe and Ewald, 2019; Marinazzo et al., 2019).

We demonstrated with several quality measures (Grova et al., 2006b; Haufe et al., 2008; Paz-Linares et al., 2017; Van de Steen et al., 2019; Wang et al., 2019) that cESI estimator distortions in actual data are more larger than expected for some state-of-the-art methods, such as sELORETA and LCMV. Our Simulated-EEG vs. MEG and EEG vs. ECoG validation method benchmarked the cESI distortions in the 10–20 EEG system (19 channels), which is considered the lower bound for all ESI. Therefore, we can infer that the behavior of our methods in EEG or other techniques that provide denser recordings can only improve.

The human Simulated-EEG vs. MEG validation method produces measurements of the cESI estimator distortions very similar to those of the EEG vs. ECoG simultaneously recorded in the macaque (Nagasaka et al., 2011; Wang et al., 2019). This similarity strongly supports the possibility of effectively measuring the distortions expected in many real-life scenarios with our Simulated-EEG vs. MEG design, which is easily replicable for other MEG acquisitions or different experimental conditions (not only resting state).

Our results indicate that ssSBL produces less cESI distortions than sELORETA and sLCMV, according to all leakage measures computed for the Simulated-EEG vs. MEG. These results were supported by cESI obtained from low-density EEG recordings in macaque compared to more fine-grained cESI obtained from simultaneously recorded high-density ECoG.

Conclusion

In this manuscript, we leveraged Bayesian theory to investigate the benefits and shortcomings of facing cESI with a large family of methods. We show that cESI must be taken care of in severely ill-conditioned and high-dimensional settings such as the ones dealt with here. In such settings, achieving cESI via exact Bayesian MAP2 methods is unfeasible and requires approximations. We have introduced quasilinearity, a reasonable cESI assumption, leading to the joint MAP that is feasible via the variational approximation to MAP. Our implementation, ssSBL specifies a priori to diminish CST distortions and exhibits good properties, outperforming state-of-the-art methods. The Bayesian theory and methods presented here could potentially be applied to signal processing and imaging other biological phenomena described by the cross-spectrum.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author/s.

Author contributions

DP-L and EG-M contributed in developing the theory, mathematical demonstrations, design of methods, implementation, and paper preparation. AA-G, ML, and YW prepared part of the data, code, and results used in this work. MV-H, QW, and EM-M prepared the macaque data for experimental confirmation and contributed to the theory and methods. MB-V, JB-B, and MV-S as senior authors analyzed the results. PV-S as senior author, introduced the theoretical background, CNS program of and guided this work. All authors contributed to the article and approved the submitted version.

Funding

This work was supported in part by the Chengdu MOST grant of 2022 under funding No. GH02-00042-HZ and the CNS program of the University of Electronic Sciences and Technology of China (UESTC) under funding No. Y0301902610100201.

Acknowledgments

We thank Prof. F. Xavier Castellanos for his assistance in reviewing and editing this manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2023.978527/full#supplementary-material

References

Andersen, E. B. (1970). Sufficiency and exponential families for discrete sample spaces. J. Am. Stat. Assoc. 65, 1248–1255. doi: 10.1080/01621459.1970.10481160

PubMed Abstract | CrossRef Full Text | Google Scholar

Andrews, D. F. (1974). Scale mixtures of normal distributions. J. R. Stat. Soc. Ser. B 148, 148–162.

Google Scholar

Auranen, T., Nummenmaa, A., Hämäläinen, M. S., Jääskeläinen, I. P., Lampinen, J., Vehtari, A., et al. (2005). Bayesian analysis of the neuromagnetic inverse problem with l(p)-norm priors. Neuroimage 26, 870–884. doi: 10.1016/j.neuroimage.2005.02.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Babacan, S. D., Mancera, L., Molina, R., and Katsaggelos, A. K. (2009). Bayesian compressive sensing using non-convex priors. Eur. Signal Process. Conf. 19, 53–63. doi: 10.1109/TIP.2009.2032894

PubMed Abstract | CrossRef Full Text | Google Scholar

Baillet, S., Mosher, J. C., and Leahy, R. M. (2001). Electromagnetic brain mapping. IEEE Signal Process. Mag. 18, 14–30. doi: 10.1109/79.962275

CrossRef Full Text

Blei, D. M., and Jordan, M. I. (2006). Variational inference for dirichlet process mixtures. Bayesian Anal. 1, 121–143. doi: 10.1214/06-BA104

CrossRef Full Text | Google Scholar

Böhning, D., and Seidel, W. (2003). Recent developments in mixture models. Comput. Stat. Data Anal. 41, 349–357. doi: 10.1016/S0167-9473(02)00161-5

CrossRef Full Text | Google Scholar

Bourbaki, N. (2013). Topological Vector Spaces: Chapters 1–5. Berlin: Springer Science & Business Media.

Google Scholar

Bradley, A., Yao, J., Dewald, J., and Richter, C.-P. (2016). Evaluation of electroencephalography source localization algorithms with multiple cortical sources. PLoS ONE 11, e0147266. doi: 10.1371/journal.pone.0147266

PubMed Abstract | CrossRef Full Text | Google Scholar

Brillinger, D. R. (1965). An introduction to polyspectra. Ann. Math. Stat. 1351–1374. doi: 10.1214/aoms/1177699896

CrossRef Full Text | Google Scholar

Brillinger, D. R. (2001). Time Series: Data Analysis and Theory. Philadelphia, PA: SIAM.

PubMed Abstract | Google Scholar

Brillinger, D. R. (2012). “Asymptotic normality of finite fourier transforms of stationary generalized processes,” in Selected Works of David Brillinger (Springer), 65–72. doi: 10.1007/978-1-4614-1344-8_6

CrossRef Full Text | Google Scholar

Brillinger, D. R., and Rosenblatt, M. (1967). Asymptotic theory of estimates of kth-order spectra. Proc. Nat. Acad. Sci. U. S. A. 57, 206–210. doi: 10.1073/pnas.57.2.206

PubMed Abstract | CrossRef Full Text | Google Scholar

Brookes, M. J., Hale, J. R., Zumer, J. M., Stevenson, C. M., Francis, S. T., Barnes, G. R., et al. (2011a). Measuring functional connectivity using meg: methodology and comparison with FcMRI. Neuroimage 56, 1082–1104. doi: 10.1016/j.neuroimage.2011.02.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Brookes, M. J., Woolrich, M., Luckhoo, H., Price, D., Hale, J. R., Stephenson, M. C., et al. (2011b). Investigating the electrophysiological basis of resting state networks using magnetoencephalography. Proc. Natl. Acad. Sci. U. S. A. 108, 16783–16788. doi: 10.1073/pnas.1112685108

PubMed Abstract | CrossRef Full Text | Google Scholar

Brookes, M. J., Woolrich, M. W., and Barnes, G. R. (2012). Measuring functional connectivity in meg: a multivariate approach insensitive to linear source leakage. Neuroimage 63, 910–912. doi: 10.1016/j.neuroimage.2012.03.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruns, A. (2004). Fourier-, Hilbert- and Wavelet-based signal analysis: are they really different approaches? J. Neurosci. Methods 137, 237. doi: 10.1016/j.jneumeth.2004.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Bryant, M., and Sudderth, E. (2012). Truly nonparametric online variational inference for hierarchical dirichlet processes. Adv. Neural Inf. Process. Syst. 25, 2699–2707. doi: 10.5555/2999325.2999436

CrossRef Full Text | Google Scholar

Burle, B., Spieser, L., Roger, C., Casini, L., Hasbroucq, T., and Vidal, F. (2015). Spatial and temporal resolutions of EEG: is it really black and white? A scalp current density view. Int. J. Psychophysiol. 97, 210–220. doi: 10.1016/j.ijpsycho.2015.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsáki, G., Logothetis, N., and Singer, W. (2013). Scaling brain size, keeping timing: evolutionary preservation of brain rhythms. Neuron 80, 751–764. doi: 10.1016/j.neuron.2013.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Candes, E., and Recht, B. (2012). Exact matrix completion via convex optimization. Commun. ACM 55, 111–119. doi: 10.1145/2184319.2184343

CrossRef Full Text | Google Scholar

Casella, G., and Berger, R. L. (2021). Statistical Inference. Boston, MA: Cengage Learning.

Google Scholar

Casella, G., Ghosh, M., Gill, J., and Kyung, M. (2010). Penalized regression, standard errors, and bayesian lassos. Bayesian Anal. 5, 369–411. doi: 10.1214/10-BA607

CrossRef Full Text | Google Scholar

Chen, K., Dong, H., and Chan, K.-S. (2013). Reduced rank regression via adaptive nuclear norm penalization. Biometrika 100, 901–920. doi: 10.1093/biomet/ast036

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Li, Y., Yu, W., Luo, L., Chen, W., and Toumoulin, C. (2011). Joint-MAP tomographic reconstruction with patch similarity based mixture prior model. Multiscale Model. Simul. 9, 1399–1419. doi: 10.1137/100814184

CrossRef Full Text | Google Scholar

Clark, V. P., Fan, S., and Hillyard, S. A. (1994). Identification of early visual evoked potential generators by retinotopic and topographic analyses. Hum. Brain Mapp. 2, 170–187. doi: 10.1002/hbm.460020306

CrossRef Full Text | Google Scholar

Colclough, G. L., Brookes, M. J., Smith, S. M., and Woolrich, M. W. (2015). A symmetric multivariate leakage correction for MEG connectomes. Neuroimage 117, 439–448. doi: 10.1016/j.neuroimage.2015.03.071

PubMed Abstract | CrossRef Full Text | Google Scholar

Csilléry, K., Blum, M. G. B., Gaggiotti, O. E., and François, O. (2010). Approximate Bayesian computation (ABC) in practice. Trends Ecol. Evol. 25, 410–418. doi: 10.1016/j.tree.2010.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Daunizeau, J., David, O., and Stephan, K. E. (2011). Dynamic causal modelling: a critical review of the biophysical and statistical foundations. Neuroimage 58, 312–322. doi: 10.1016/j.neuroimage.2009.11.062

PubMed Abstract | CrossRef Full Text | Google Scholar

Daunizeau, J., and Friston, K. J. (2007). A mesostate-space model for EEG and MEG. Neuroimage 38, 67–81. doi: 10.1016/j.neuroimage.2007.06.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Daunizeau, J., Laufs, H., and Friston, K. J. (2010). EEG-FMRI information fusion: biophysics and data analysis. EEGFMRI Physiol. Basis Techn. Applic. 511–26. doi: 10.1007/978-3-540-87919-0_25

CrossRef Full Text | Google Scholar

Davis, L. M., Collings, I. B., and Hoeher, P. (2001). Joint MAP equalization and channel estimation for frequency-selective and frequency-flat fast-fading channels. IEEE Trans. Commun. 49, 2106–2114. doi: 10.1109/26.974257

CrossRef Full Text | Google Scholar

Deco, G., Jirsa, V. K., Robinson, P. A., Breakspear, M., and Friston, K. (2008). The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Comput. Biol. 4, e1000092. doi: 10.1371/journal.pcbi.1000092

PubMed Abstract | CrossRef Full Text | Google Scholar

Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). maximum likelihood from incomplete data via the EM algorithm. J. R. Stat. Soc. Ser. B 39, 1–22. doi: 10.1111/j.2517-6161.1977.tb01600.x

CrossRef Full Text | Google Scholar

Ding, C., Zhou, D., He, X., and Zha, H. (2006). “R 1-Pca: rotational invariant l 1-norm principal component analysis for robust subspace factorization,” in Proceedings of the 23rd International Conference on Machine Learning (Pittsburgh, PA), 281–88. doi: 10.1145/1143844.1143880

CrossRef Full Text | Google Scholar

Dunford, N., and Schwartz, J. T. (1988). Linear Operators, Part 1: General Theory. Vol. 10. Hoboken, NJ: John Wiley & Sons.

Google Scholar

Duvenaud, D., Maclaurin, D., and Adams, R. (2016). “Early stopping as nonparametric variational inference,” in Artificial Intelligence and Statistics (PMLR) , 1070–1077.

Google Scholar

Eichele, T., Specht, K., Moosmann, M., Jongsma, M. L. A., Quiroga, R. Q., Nordby, H., et al. (2005). Assessing the spatiotemporal evolution of neuronal activation with single-trial event-related potentials and functional MRI. Proc. Nat. Acad. Sci. U. S. A. 102, 17798–17803. doi: 10.1073/pnas.0505508102

PubMed Abstract | CrossRef Full Text | Google Scholar

Engel, A. K., Fries, P., and Singer, W. (2001). Dynamic predictions: oscillations and synchrony in top–down processing. Nat. Rev. Neurosci. 2, 704–716. doi: 10.1038/35094565

PubMed Abstract | CrossRef Full Text | Google Scholar

Eyink, G. L., Restrepo, J. M., and Alexander, F. J. (2004). A mean field approximation in data assimilation for nonlinear dynamics. Phys. D Nonlinear Phenomena 195, 347–368. doi: 10.1016/j.physd.2004.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Faes, L., Erla, S., and Nollo, G. (2012). Measuring connectivity in linear multivariate processes: definitions, interpretation, and practical analysis. Comput. Math. Methods Med. 2012, 140513. doi: 10.1155/2012/140513

PubMed Abstract | CrossRef Full Text | Google Scholar

Faes, L., and Nollo, G. (2011). Multivariate frequency domain analysis of causal interactions in physiological time series. Biomed. Eng. Trends Electron. Commun. Softw. 8, 403–428. doi: 10.5772/13065

CrossRef Full Text | Google Scholar

Faes, L., Stramaglia, S., and Marinazzo, D. (2017). On the interpretability and computational reliability of frequency-domain Granger causality. F1000 Res. 6, 1710. doi: 10.12688/f1000research.12694.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, K. (1951). Maximum properties and inequalities for the eigenvalues of completely continuous operators. Proc. Nat. Acad. Sci. U. S. A. 37, 760–766. doi: 10.1073/pnas.37.11.760

PubMed Abstract | CrossRef Full Text | Google Scholar

Frauscher, B., Ellenrieder, N. V., Zelmann, R., DoleŽalová, I., Minotti, L., Olivier, A., et al. (2018). Atlas of the normal intracranial electroencephalogram: neurophysiological awake activity in different cortical areas. Brain 141, 1130–1144. doi: 10.1093/brain/awy035

PubMed Abstract | CrossRef Full Text | Google Scholar

Freeman, W. J. (1975). Mass Action in the Nervous System : Examination of the Neurophysiological Basis of Adaptive Behavior through the EEG. New York, NY; London: Academic Press. Available online at: https://librarysearch.cardiff.ac.uk/openurl/44WHELF_CAR/44WHELF_CAR:44WHELF_CAR_VU1??u.ignore_date_coverage=true&rft.mms_id=9911670637102420 (accessed February 01, 2023).

Google Scholar

Friston, K., Harrison, L., Daunizeau, J., Kiebel, S., Phillips, C., Trujillo-Barreto, N., et al. (2008). Multiple sparse priors for the M/EEG inverse problem. Neuroimage 39, 1104–1120. doi: 10.1016/j.neuroimage.2007.09.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K., Henson, R., Phillips, C., and Mattout, J. (2006). Bayesian estimation of evoked and induced responses. Hum. Brain Mapp. 27, 722–735. doi: 10.1002/hbm.20214

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K., Mattout, J., Trujillo-Barreto, N., Ashburner, J., and Penny, W. (2007). Variational free energy and the laplace approximation. Neuroimage 34, 220–234. doi: 10.1016/j.neuroimage.2006.08.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J. (2009). Modalities, modes, and models in functional neuroimaging. Science 326, 399–403. doi: 10.1126/science.1174521

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J., Bastos, A., Litvak, V., Stephan, K. E., Fries, P., and Moran, R. J. (2012). DCM for complex-valued data: cross-spectra, coherence and phase-delays. Neuroimage 59, 439–455. doi: 10.1016/j.neuroimage.2011.07.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Friston, K. J., Trujillo-Barreto, N., and Daunizeau, J. (2008). DEM: a variational treatment of dynamic systems. Neuroimage 41, 849–885. doi: 10.1016/j.neuroimage.2008.02.054

PubMed Abstract | CrossRef Full Text | Google Scholar

Gershman, S., Hoffman, M., and Blei, D. (2012). Nonparametric variational inference. ArXiv Preprint ArXiv:1206.4665. 53, 1–12.

Google Scholar

Gershman, S. J., and Blei, D. M. (2012). A tutorial on Bayesian nonparametric models. J. Math. Psychol. 56, 1–12. doi: 10.1016/j.jmp.2011.08.004

CrossRef Full Text | Google Scholar

Ghahramani, Z., and Beal, M. J. (2001). Propagation algorithms for variational bayesian learning. Adv. Neural Inform. Process. Syst. 3, 507–513.

Google Scholar

Ghahramani, Z., and Beal, M. M. J. (2000). Variational inference for Bayesian mixtures of factor analysers. Adv. Neural Inform. Process. Syst. 12, 449–455.

Google Scholar

Golub, G. H., and Van Loan, C. F. (2013). Matrix Computations. Charles Village, MA: JHU Press. doi: 10.56021/9781421407944

CrossRef Full Text | Google Scholar

Gramfort, A., Kowalski, M., and Hämäläinen, M. (2012). Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods. Phys. Med. Biol. 57, 1937. doi: 10.1088/0031-9155/57/7/1937

PubMed Abstract | CrossRef Full Text | Google Scholar

Gramfort, A., Strohmeier, D., Haueisen, J., Hämäläinen, M. S., and Kowalski, M. (2013). Time-frequency mixed-norm estimates: sparse M/EEG imaging with non-stationary source activations. Neuroimage 70, 410–422. doi: 10.1016/j.neuroimage.2012.12.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Grave de Peralta Menendez, R., Murray, M. M., Michel, C. M., Martuzzi, R., and Andino, S. L. G. (2004). Electrical neuroimaging based on biophysical constraints. Neuroimage 21, 527–539. doi: 10.1016/j.neuroimage.2003.09.051

PubMed Abstract | CrossRef Full Text

Grech, R., Cassar, T., Muscat, J., Camilleri, K. P., Fabri, S. G., Zervakis, M., et al. (2008). Review on solving the inverse problem in EEG source analysis. J. Neuroeng. Rehabil. 5, 1–33. doi: 10.1186/1743-0003-5-25

PubMed Abstract | CrossRef Full Text | Google Scholar

Grova, C., Daunizeau, J., Kobayashi, E., Bagshaw, A. P., Lina, J. M., Dubeau, F., et al. (2008). Concordance between distributed EEG source localization and simultaneous EEG-FMRI studies of epileptic spikes. Neuroimage 39, 755–774. doi: 10.1016/j.neuroimage.2007.08.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Grova, C., Daunizeau, J., Lina, J. M., Bénar, C. G., Benali, H., and Gotman, J. (2006a). Evaluation of EEG localization methods using realistic simulations of interictal spikes. Neuroimage 29, 734–753. doi: 10.1016/j.neuroimage.2005.08.053

PubMed Abstract | CrossRef Full Text | Google Scholar

Grova, C., Makni, S., Flandin, G., Ciuciu, P., Gotman, J., and Poline, J. B. (2006b). Anatomically informed interpolation of FMRI data on the cortical surface. Neuroimage 31, 1475–1486. doi: 10.1016/j.neuroimage.2006.02.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Hadamard, J., and Morse, P. M. (1953). Lectures on Cauchy's problem in linear partial differential equations. Phys. Today 6, 18–18. doi: 10.1063/1.3061337

CrossRef Full Text | Google Scholar

Hallez, H., Vanrumste, B., Grech, R., Muscat, J., Clercq, W. D., Vergult, A., et al. (2007). Review on solving the forward problem in EEG source analysis. J. Neuroeng. Rehabil. 4, 46. doi: 10.1186/1743-0003-4-46

PubMed Abstract | CrossRef Full Text | Google Scholar

Hämäläinen, M., Hari, R., Ilmoniemi, R. J., Knuutila, J., and Lounasmaa, O. V. (1993). Magnetoencephalography—theory, instrumentation, and applications to noninvasive studies of the working human brain. Rev. Mod. Phys. 65, 413. doi: 10.1103/RevModPhys.65.413

CrossRef Full Text | Google Scholar

Hämäläinen, M. S., and Ilmoniemi, R. J. (1994). Interpreting magnetic fields of the brain: minimum norm estimates. Med. Biol. Eng. Comput. 32, 35–42. doi: 10.1007/BF02512476

PubMed Abstract | CrossRef Full Text | Google Scholar

Hancock, G. R., and Samuelsen, K. M. (2007). Advances in Latent Variable Mixture Models. Charlotte, NC: IAP.

Google Scholar

Harrison, L. M., Penny, W., Daunizeau, J., and Friston, K. J. (2008). Diffusion-based spatial priors for functional magnetic resonance images. Neuroimage 41, 408–423. doi: 10.1016/j.neuroimage.2008.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Haufe, S., and Ewald, A. (2019). A simulation framework for benchmarking EEG-based brain connectivity estimation methodologies. Brain Topogr. 32, 625–642. doi: 10.1007/s10548-016-0498-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Haufe, S., Nikulin, V. V., Müller, K. R., and Nolte, G. (2013). A critical assessment of connectivity measures for EEG data: a simulation study. Neuroimage 64, 120–133. doi: 10.1016/j.neuroimage.2012.09.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Haufe, S., Nikulin, V. V., Ziehe, A., Müller, K. R., and Nolte, G. (2008). Combining sparsity and rotational invariance in EEG/MEG source reconstruction. Neuroimage 42, 726–738. doi: 10.1016/j.neuroimage.2008.04.246

PubMed Abstract | CrossRef Full Text | Google Scholar

Hauk, O. (2004). Keep it simple: a case for using classical minimum norm estimation in the analysis of EEG and MEG data. Neuroimage 21, 1612–1621. doi: 10.1016/j.neuroimage.2003.12.018

PubMed Abstract | CrossRef Full Text | Google Scholar

He, B., Astolfi, L., Valdes-Sosa, P. A., Marinazzo, D., Palva, S. O., Benar, C.-G., et al. (2019). Electrophysiological brain connectivity: theory and implementation. IEEE Trans. Biomed. Eng. 66, 2115–2137. doi: 10.1109/TBME.2019.2913928

PubMed Abstract | CrossRef Full Text | Google Scholar

He, B., Sohrabpour, A., Brown, E., and Liu, Z. (2018). Electrophysiological source imaging: a noninvasive window to brain dynamics. Annu. Rev. Biomed. Eng. 20:171–196. doi: 10.1146/annurev-bioeng-062117-120853

PubMed Abstract | CrossRef Full Text | Google Scholar

Henson, R. N., Wakeman, D. G., Litvak, V., and Friston, K. J. (2011). A parametric empirical bayesian framework for the EEG/MEG inverse problem: generative models for multi-subject and multi-modal integration. Front. Hum. Neurosci. 5, 76. doi: 10.3389/fnhum.2011.00076

PubMed Abstract | CrossRef Full Text | Google Scholar

Hindriks, R. (2020). A methodological framework for inverse-modeling of propagating cortical activity using MEG/EEG. Neuroimage 223:117345. doi: 10.1016/j.neuroimage.2020.117345

PubMed Abstract | CrossRef Full Text | Google Scholar

Hipp, J. F., Hawellek, D. J., Corbetta, M., Siegel, M., and Engel, A. K. (2012). Large-scale cortical correlation structure of spontaneous oscillatory activity. Nat. Neurosci. 15, 884–890. doi: 10.1038/nn.3101

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoerl, A. E., and Kennard, R. W. (1970). Ridge regression: biased estimation for nonorthogonal problems. Technometrics 12, 55–67. doi: 10.1080/00401706.1970.10488634

CrossRef Full Text | Google Scholar

Hsiao, T., Rangarajan, A., and Gindi, G. (1998). “Joint-MAP reconstruction/segmentation for transmission tomography using mixture-models as priors,” in 1998 IEEE Nuclear Science Symposium Conference Record 1998. IEEE Nuclear Science Symposium and Medical Imaging Conference (Cat. No. 98CH36255), Vol. 3 (Toronto, ON: IEEE), 1689–1693.

Google Scholar

Hsiao, T., Rangarajan, A., and Gindi, G. (2002). Joint-MAP Bayesian tomographic reconstruction with a gamma-mixture prior. IEEE Trans. Image Process. 11, 1466–1477. doi: 10.1109/TIP.2002.806254

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, Y., Zhang, D., Ye, J., Li, X., and He, X. (2012). Fast and accurate matrix completion via truncated nuclear norm regularization. IEEE Trans. Pattern Anal. Mach. Intell. 35, 2117–2130. doi: 10.1109/TPAMI.2012.271

PubMed Abstract | CrossRef Full Text | Google Scholar

Jirsa, V., and Haken, H. (1996). Field theory of electromagnetic brain activity. Phys. Rev. Lett. 77, 960–963. doi: 10.1103/PhysRevLett.77.960

PubMed Abstract | CrossRef Full Text | Google Scholar

Jirsa, V. K. (2004). Connectivity and dynamics of neural information processing. Neuroinformatics 2, 183–204. doi: 10.1385/NI:2:2:183

PubMed Abstract | CrossRef Full Text | Google Scholar

Jirsa, V. K., and Haken, H. (1997). A derivation of a macroscopic field theory of the brain from the quasi-microscopic neural dynamics. Phys. D Nonlinear Phenomena 99, 503–526. doi: 10.1016/S0167-2789(96)00166-2

CrossRef Full Text | Google Scholar

Kadanoff, L. P. (2009). More is the same; phase transitions and mean field theories. J. Stat. Phys. 137, 777–797. doi: 10.1007/s10955-009-9814-1

CrossRef Full Text | Google Scholar

Kaplan, A., and Tichatschke, R. (1998). Proximal point methods and nonconvex optimization. J. Glob. Optimiz. 13, 389–406. doi: 10.1023/A:1008321423879

CrossRef Full Text | Google Scholar

Kim, E., Lee, M., and Oh, S. (2015). “Elastic-net regularization of singular values for robust subspace learning,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (Boston, MA), 915–923. doi: 10.1109/CVPR.2015.7298693

PubMed Abstract | CrossRef Full Text | Google Scholar

Kindermann, R., Snell, J. L., Society, A. M., and Collection, K. M. R. (1980). Markov Random Fields and Their Applications. Contemporary Mathematics - American Mathematical Society. American Mathematical Society. Available online at: https://books.google.co.uk/books?id=NeVQAAAAMAAJ (accessed February 01, 2023). doi: 10.1090/conm/001

CrossRef Full Text

Knösche, T. R., and Haueisen, J. (2022). EEG/MEG Source Reconstruction: Textbook for Electro-and Magnetoencephalography. Midtown Manhattan, NY: Springer. doi: 10.1007/978-3-030-74918-7

CrossRef Full Text | Google Scholar

Kobayashi, K., Yoshinaga, H., Oka, M., Ohtsuka, Y., and Gotman, J. (2003). A simulation study of the error in dipole source localization for EEG spikes with a realistic head model. Clin. Neurophysiol. 114, 1069–1078. doi: 10.1016/S1388-2457(03)00064-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Kowalski, M., and Torrésani, B. (2009a). “Structured sparsity: from mixed norms to structured shrinkage,” in SPARS'09-Signal Processing with Adaptive Sparse Structured Representations (Saint-Malo).

Google Scholar

Kowalski, M., and Torrésani, B. (2009b). Sparsity and persistence: mixed norms provide simple signal models with dependent coefficients. Signal Image Video Process. 3, 251–264. doi: 10.1007/s11760-008-0076-1

CrossRef Full Text | Google Scholar

Lafferty, J., McCallum, A., and Pereira, F. C. N. (2001). “Conditional random fields: Probabilistic models for segmenting and labeling sequence data,” in ICML '01: Proceedings of the Eighteenth International Conference on Machine Learning (San Francisco, CA), 282–289. doi: 10.5555/645530.655813

CrossRef Full Text | Google Scholar

Landau, L. D., and Lifshitz, E. M. (1980). “Chapter III - the gibbs distribution,” in Statistical Physics, 3rd Edn., eds L. D. Landau and E. M. Lifshitz (Oxford: Butterworth-Heinemann), 79–110. doi: 10.1016/B978-0-08-057046-4.50010-5

CrossRef Full Text

Larson-Prior, L. J., Oostenveld, R., Penna, S. D., Michalareas, G., Prior, F., Babajani-Feremi, A., et al. (2013). Adding dynamics to the human connectome project with MEG. Neuroimage 80:190–201. doi: 10.1016/j.neuroimage.2013.05.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Boudec, J.-Y., McDonald, D., and Mundinger, J. (2007). “A generic mean field convergence result for systems of interacting objects,” in Fourth International Conference on the Quantitative Evaluation of Systems (QEST 2007) (Edinburgh: IEEE), 3–18. doi: 10.1109/QEST.2007.8

CrossRef Full Text | Google Scholar

Le Van Quyen, M., and Bragin, A. (2007). Analysis of dynamic brain oscillations: methodological advances. Trends Neurosci. 30, 365–373. doi: 10.1016/j.tins.2007.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Lei, X., Ostwald, D., Hu, J., Qiu, C., Porcaro, C., Bagshaw, A. P., et al. (2011). Multimodal functional network connectivity: an EEG-FMRI fusion in network space. PLoS ONE 6, e24642. doi: 10.1371/journal.pone.0024642

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Li, R., and Wu, R. (2011). Bayesian group LASSO for nonparametric varying-coefficient. Group 16802, 1–19.

PubMed Abstract | Google Scholar

Li, Q., and Lin, N. (2010). The Bayesian elastic net. Bayesian Anal. 5, 151–170. doi: 10.1214/10-BA506

CrossRef Full Text | Google Scholar

Li, Z., and Tian, T. S. (2011). A spatio-temporal solution for the EEG/MEG inverse problem using group penalization methods. Stat. Interface 4, 521–533. doi: 10.4310/SII.2011.v4.n4.a10

CrossRef Full Text | Google Scholar

Lindsay, B. G. (1995). Mixture models: theory, geometry, and applications. NSF-CBMS Regional Conf. Ser. Probab. Statist. 5, 163. doi: 10.1214/cbms/1462106013

CrossRef Full Text | Google Scholar

Liu, C., and Rubin, D. B. (1994). The ECME algorithm: a simple extension of EM and ECM with faster monotone convergence. Biometrika 81, 633–648. doi: 10.1093/biomet/81.4.633

CrossRef Full Text | Google Scholar

Lopes da Silva, F. (2013). EEG and MEG: relevance to neuroscience. Neuron 80, 1112–1128. doi: 10.1016/j.neuron.2013.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes da Silva, F. H., Hoeks, A., Smits, H., and Zetterberg, L. H. (1974). Model of brain rhythmic activity - the alpha-rhythm of the thalamus. Kybernetik 15, 27–37. doi: 10.1007/BF00270757

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes da Silva, F. H., Wieringa, H. J., and Peters, M. J. (1991). Source localization of EEG versus MEG: empirical comparison using visually evoked responses and theoretical considerations. Brain Topogr. 4, 133–142. doi: 10.1007/BF01132770

PubMed Abstract | CrossRef Full Text | Google Scholar

MacKay, D. J. C. (1999). Comparison of approximate methods for handling hyperparameters. Neural Comput. 11, 1035–1068. doi: 10.1162/089976699300016331

CrossRef Full Text | Google Scholar

MacKay, D. J. C. (2003). Information Theory, Inference and Learning Algorithms. Cambridge: Cambridge University Press.

Google Scholar

MacKay, D. J. C. J. C. (1998). Choice of basis for laplace approximation. Mach. Learn. 33, 77–86. doi: 10.1023/A:1007558615313

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahjoory, K., Nikulin, V. V., Botrel, L., Linkenkaer-Hansen, K., Fato, M. M., and Haufe, S. (2017). Consistency of EEG source localization and connectivity estimates. Neuroimage 152:590–601. doi: 10.1016/j.neuroimage.2017.02.076

PubMed Abstract | CrossRef Full Text | Google Scholar

Makeig, S. (2002). Response: event-related brain dynamics - unifying brain electrophysiology. Trends Neurosci. 25, 390. doi: 10.1016/S0166-2236(02)02198-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Makeig, S., Delorme, A., Westerfield, M., Jung, T. P., Townsend, J., Courchesne, E., et al. (2004). Electroencephalographic brain dynamics following manually responded visual targets. PLoS Biol. 2, e176. doi: 10.1371/journal.pbio.0020176

PubMed Abstract | CrossRef Full Text | Google Scholar

Makeig, S., Westerfield, M., Jung, T. P., Covington, J., Townsend, J., Sejnowski, T. J., et al. (1999). Functionally independent components of the late positive event-related potential during visual spatial attention. J. Neurosci. 19, 2665–2680. doi: 10.1523/JNEUROSCI.19-07-02665.1999

PubMed Abstract | CrossRef Full Text | Google Scholar

Mantini, D., Perrucci, M. G., Gratta, C. D., Romani, G. L., and Corbetta, M. (2007). Electrophysiological signatures of resting state networks in the human brain. Proc. Nat. Acad. Sci. U. S. A. 104, 13170–13175. doi: 10.1073/pnas.0700668104

PubMed Abstract | CrossRef Full Text | Google Scholar

Marinazzo, D., Riera, J. J., Marzetti, L., Astolfi, L., Yao, D., and Sosa, P. A. V. (2019). Controversies in EEG source imaging and connectivity: modeling, validation, benchmarking. Brain Topogr. 32, 527–529. doi: 10.1007/s10548-019-00709-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Marzetti, L., Gratta, C. D., and Nolte, G. (2008). Understanding brain connectivity from EEG data by identifying systems composed of interacting sources. Neuroimage 42, 87–98. doi: 10.1016/j.neuroimage.2008.04.250

PubMed Abstract | CrossRef Full Text | Google Scholar

Mattout, J., Phillips, C., Penny, W. D., Rugg, M. D., and Friston, K. J. (2006). MEG source localization under multiple constraints: an extended bayesian framework. Neuroimage 30, 753–767. doi: 10.1016/j.neuroimage.2005.10.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Maurer, K., and Dierks, T. (2012). Atlas of Brain Mapping: Topographic Mapping of EEG and Evoked Potentials. Berlin; Heidelberg: Springer Science & Business Media.

Google Scholar

McLachlan, G. J., and Basford, K. E. (1988). Mixture Models: Inference and Applications to Clustering. Vol. 38. New York, NY: M. Dekker New York. doi: 10.2307/2348072

CrossRef Full Text | Google Scholar

Moran, R., Pinotsis, D. a., and Friston, K. (2013). Neural masses and fields in dynamic causal modeling. Front. Comput. Neurosci. 7, 57. doi: 10.3389/fncom.2013.00057

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagasaka, Y., Shimoda, K., and Fujii, N. (2011). Multidimensional recording (MDR) and data sharing: an ecological open research and educational platform for neuroscience. PLoS ONE 6, e22561. doi: 10.1371/journal.pone.0022561

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, T., and Bonilla, E. (2013). “Efficient variational inference for gaussian process regression networks,” in Artificial Intelligence and Statistics (Scottsdale, AZ: PMLR), 472–480.

Google Scholar

Niedermeyer, E., and da Silva, F. H. L. (2005). Electroencephalography: Basic Principles, Clinical Applications, and Related Fields. LWW Doody's All Reviewed Collection. Lippincott Williams & Wilkins. Available online at: https://books.google.com.sg/books?id=tndqYGPHQdEC (accessed February 01, 2023).

Google Scholar

Nolte, G., Galindo-Leon, E., Li, Z., Liu, X., and Engel, A. K. (2020). Mathematical relations between measures of brain connectivity estimated from electrophysiological recordings for gaussian distributed data. Front. Neurosci. 14, 577574. doi: 10.3389/fnins.2020.577574

PubMed Abstract | CrossRef Full Text | Google Scholar

Nummenmaa, A., Auranen, T., Hämäläinen, M. S., Jääskeläinen, I. P., Lampinen, J., Sams, M., et al. (2007). Hierarchical Bayesian estimates of distributed MEG sources: theoretical aspects and comparison of variational and MCMC methods. Neuroimage 35, 669–685. doi: 10.1016/j.neuroimage.2006.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunez, P. L. (1974). The brain wave equation: a model for the EEG. Math. Biosci. 21, 279–97. doi: 10.1016/0025-5564(74)90020-0

CrossRef Full Text | Google Scholar

Nunez, P. L., Silberstein, R. B., Cadusch, P. J., Wijesinghe, R. S., Westdorp, A. F., and Srinivasan, R. (1994). A theoretical and experimental study of high resolution eeg based on surface laplacians and cortical imaging. Electroencephalogr. Clin. Neurophysiol. 90, 40–57. doi: 10.1016/0013-4694(94)90112-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Nunez, P. L., and Srinivasan, R. (2006). Electric Fields of the Brain: The Neurophysics of EEG, 2nd Edn. New York, NY: Oxford University Press. doi: 10.1093/acprof:oso/9780195050387.001.0001

PubMed Abstract | CrossRef Full Text | Google Scholar

Palva, J. M., Wang, S. H., Palva, S., Zhigalov, A., Monto, S., Brookes, M. J., et al. (2018). Ghost interactions in MEG/EEG source space: a note of caution on inter-areal coupling measures. Neuroimage 173, 632–643. doi: 10.1016/j.neuroimage.2018.02.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Papadopoulou, M., Friston, K., and Marinazzo, D. (2019). Estimating directed connectivity from cortical recordings and reconstructed sources. Brain Topogr. 32, 741–752. doi: 10.1007/s10548-015-0450-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, T., and Casella, G. (2008). The Bayesian lasso. J. Am. Stat. Assoc. 103, 681–686. doi: 10.1198/016214508000000337

CrossRef Full Text | Google Scholar

Pascual-Marqui, R. D., Pascual-Montano, A. D., Lehmann, D., Kochi, K., Esslen, M., Jancke, L., et al. (2006). Exact low resolution brain electromagnetic tomography (ELORETA). Neuroimage 31 (Suppl. 1):S86. doi: 10.1016/S1053-8119(08)70001-6

CrossRef Full Text | Google Scholar

Paz-Linares, D., Gonzalez-Moreira, E., Areces-Gonzalez, A., Wang, Y., Li, M., Martinez-Montes, E., et al. (2018). Identification of oscillatory brain networks with hidden gaussian graphical spectral models of EEG/MEG. ArXiv [Preprint]. arXiv:1810.01174. doi: 10.48550/arXiv.1810.01174

CrossRef Full Text | Google Scholar

Paz-Linares, D., Vega-Hernández, M., Rojas-López, P. A., Valdés-Hernández, P. A., Martínez-Montes, E., and Valdés-Sosa, P. A. (2017). Spatio temporal EEG source imaging with the hierarchical bayesian elastic net and elitist lasso models. Front. Neurosci. 11, 635. doi: 10.3389/fnins.2017.00635

PubMed Abstract | CrossRef Full Text | Google Scholar

Pearl, J. (1988). Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. New York, NY: ACM. doi: 10.1016/B978-0-08-051489-5.50008-4

CrossRef Full Text | Google Scholar

Pearl, J. (2022). “Reverend Bayes on inference engines: a distributed hierarchical approach,” in Probabilistic and Causal Inference: The Works of Judea Pearl (San Rafael, CA: Morgan & Claypool), 129–38. doi: 10.1145/3501714.3501727

CrossRef Full Text | Google Scholar

Petersen, K. B., and Pedersen, M. S. (2008). The Matrix Cookbook. Lyngby: Technical University of Denmark. Vol. 7, 510.

Google Scholar

Piastra, M. C., Nüßing, A., Vorwerk, J., Clerc, M., Engwer, C., and Wolters, C. H. (2020). A comprehensive study on electroencephalography and magnetoencephalography sensitivity to cortical and subcortical sources. Hum. Brain Mapp. 42, 978–992. doi: 10.1002/hbm.25272

PubMed Abstract | CrossRef Full Text | Google Scholar

Picton, T. W., and Hillyard, S. A. (1974). Human auditory evoked potentials. II: effects of attention. Electroencephal. Clin. Neurophysiol. 36, 191–200. doi: 10.1016/0013-4694(74)90156-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Picton, T. W., Hillyard, S. A., Krausz, H. I., and Galambos, R. (1974). Human auditory evoked potentials. I: evaluation of components. Electroencephalogr. Clin. Neurophysiol. 36:179–190. doi: 10.1016/0013-4694(74)90155-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Piotrowski, T., and Yamada, I. (2008). MV-PURE estimator: minimum-variance pseudo-unbiased reduced-rank estimator for linearly constrained ill-conditioned inverse problems. IEEE Trans. Signal Process. 56, 3408–3423. doi: 10.1109/TSP.2008.921716

CrossRef Full Text | Google Scholar

Reid, A. T., Headley, D. B., Mill, R. D., Sanchez-Romero, R., Uddin, L. Q., Marinazzo, D., et al. (2019). Advancing functional connectivity research from association to causation. Nat. Neurosci. 22, 1751–1760. doi: 10.1038/s41593-019-0510-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Riera, J. J., and Fuentes, M. E. (1998). Electric lead field for a piecewise homogeneous volume conductor model of the head. IEEE Trans. Biomed. Eng. 45, 746–753. doi: 10.1109/10.678609

PubMed Abstract | CrossRef Full Text | Google Scholar

Riesz, F. (1910). Untersuchungen über systeme integrierbarer funktionen. Math. Annalen 69, 449–497. doi: 10.1007/BF01457637

CrossRef Full Text | Google Scholar

Rosa, M. J., Daunizeau, J., and Friston, K. J. (2010). EEG-FMRI integration: a critical review of biophysical modeling and data analysis approaches. J. Integr. Neurosci. 9, 453–476. doi: 10.1142/S0219635210002512

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenblatt, M. (1956). A central limit theorem and a strong mixing condition. Proc. Natl. Acad. Sci. U. S. A. 42, 43. doi: 10.1073/pnas.42.1.43

PubMed Abstract | CrossRef Full Text | Google Scholar

Roweis, S., and Ghahramani, Z. (1999). A unifying review of linear gaussian models. Neural Comput. 11, 305–345. doi: 10.1162/089976699300016674

PubMed Abstract | CrossRef Full Text | Google Scholar

Rudin, W. (1970). Real and Complex Analysis P. 2. New York, NY: McGraw-Hill.

Google Scholar

Salmelin, R. H., and Hämäläinen, M. S. (1995). Dipole modelling of MEG rhythms in time and frequency domains. Brain Topogr. 7, 251–257. doi: 10.1007/BF01202384

PubMed Abstract | CrossRef Full Text | Google Scholar

Schatten, R. (2013). Norm Ideals of Completely Continuous Operators. Vol. 27. Midtown Manhattan, NY: Springer-Verlag.

Google Scholar

Schoffelen, J., and Gross, J. (2009). Source connectivity analysis with MEG and EEG. Hum. Brain Mapp. 30, 1857–1865. doi: 10.1002/hbm.20745

PubMed Abstract | CrossRef Full Text | Google Scholar

Stokes, P. A., and Purdon, P. L. (2017). A study of problems encountered in granger causality analysis from a neuroscience perspective. Proc. Nat. Acad. Sci. U. S. A. 201704663. doi: 10.1073/pnas.1704663114

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, T., and Zhang, C.-H. (2012). Calibrated elastic regularization in matrix completion. Adv. Neural Inf. Process. Syst. 25, 863–871.

Google Scholar

Tarantola, A. (2005). “Inverse problem theory and methods for model parameter estimation,” in Inverse Problem Theory and Methods for Model Parameter Estimation (Philadelphia, PA: SIAM). doi: 10.1137/1.9780898717921

CrossRef Full Text | Google Scholar

Tewarie, P., Liuzzi, L., O'Neill, G. C., Quinn, A. J., Griffa, A., Woolrich, M. W., et al. (2019). Tracking dynamic brain networks using high temporal resolution MEG measures of functional connectivity. Neuroimage 200, 38–50. doi: 10.1016/j.neuroimage.2019.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Tibshirani, R. (1996). Regression shrinkage and selection via the Lasso. J. R. Stat. Soc. Ser. B 58, 267–288. doi: 10.1111/j.2517-6161.1996.tb02080.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tikhonov, A. N., and Arsenin, V. I. A. (1977). Solutions of Ill-Posed Problems. Washington; New York, NY: Winston; Distributed solely by Halsted Press. Available online at: http://books.google.com/books?id=ECrvAAAAMAAJ (accessed February 01, 2023).

Google Scholar

Tipping, M. E. (2001). Sparse Bayesian learning and the relevance vector machine. J. Mach. Learn. Res. 1, 211–244. doi: 10.1162/15324430152748236

CrossRef Full Text | Google Scholar

Tirer, T., and Giryes, R. (2020). Back-projection based fidelity term for ill-posed linear inverse problems. IEEE Trans. Image Process. 29, 6164–6179. doi: 10.1109/TIP.2020.2988779

PubMed Abstract | CrossRef Full Text | Google Scholar

Valdés, P., Bosch, J., Grave, R., Hernandez, J., Riera, J., Pascual, R., et al. (1992). Frequency domain models of the EEG. Brain Topogr. 4, 309–319. doi: 10.1007/BF01135568

PubMed Abstract | CrossRef Full Text | Google Scholar

Valdes-Sosa, P. A., Sanchez-Bornot, J. M., Sotero, R. C., Iturria-Medina, Y., Aleman-Gomez, Y., Bosch-Bayard, J., et al. (2009). Model driven EEG/FMRI fusion of brain oscillations. Hum. Brain Mapp. 30, 2701–2721. doi: 10.1002/hbm.20704

PubMed Abstract | CrossRef Full Text | Google Scholar

Van de Steen, F., Faes, L., Karahan, E., Songsiri, J., Valdes-Sosa, P. A., and Marinazzo, D. (2019). Critical comments on EEG sensor space dynamical connectivity analysis. Brain Topogr. 32, 643–654. doi: 10.1007/s10548-016-0538-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Essen, D. C., Smith, S. M., Barch, D. M., Behrens, T. E. J., Yacoub, E., and Ugurbil, K. (2013). The WU-Minn human connectome project: an overview. Neuroimage 80:62–79. doi: 10.1016/j.neuroimage.2013.05.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Veen, B. D., Drongelen, W. V., Yuchtman, M., and Suzuki, A. (1997). Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans. Biomed. Eng. 44, 867–880. doi: 10.1109/10.623056

PubMed Abstract | CrossRef Full Text | Google Scholar

Varela, F., Lachaux, J.-P., Rodriguez, E., and Martinerie, J. (2001). The brainweb: phase synchronization and large-scale integration. Nat. Rev. Neurosci. 2, 229–239. doi: 10.1038/35067550

PubMed Abstract | CrossRef Full Text | Google Scholar

Vega-Hernández, M., Martínez-Montes, E., Sánchez-Bornot, J. M., Lage-Castellanos, A., and Valdés-Sosa, P. A. (2008). Penalized least squares methods for solving the EEG inverse problem. Stat. Sin. 18, 1535–1551. Available online at: https://www3.stat.sinica.edu.tw/statistica/J18N4/j18n416/j18n416.html

Google Scholar

Vidaurre, D., Abeysuriya, R., Becker, R., Quinn, A. J., Alfaro-Almagro, F., Smith, S. M., et al. (2018a). Discovering dynamic brain networks from big data in rest and task. Neuroimage 180, 646–656. doi: 10.1016/j.neuroimage.2017.06.077

PubMed Abstract | CrossRef Full Text | Google Scholar

Vidaurre, D., Hunt, L. T., Quinn, A. J., Hunt, B. A. E. E., Brookes, M. J., Nobre, A. C., et al. (2018b). Spontaneous cortical activity transiently organises into frequency specific phase-coupling networks. Nat. Commun. 9, 2987. doi: 10.1038/s41467-018-05316-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Vinck, M., and Perrenoud, Q. (2019). Layers of rhythms—from cortical anatomy to dynamics. Neuron 101, 358–360. doi: 10.1016/j.neuron.2019.01.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Valdés-Hernández, P. A., Paz-Linares, D., Bosch-Bayard, J., Oosugi, N., Komatsu, M., et al. (2019). EECoG-comp: an open source platform for concurrent EEG/ECoG comparisons—applications to connectivity studies. Brain Topogr. 32, 550–568. doi: 10.1007/s10548-019-00708-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiss, P. (1907). L'hypothèse Du Champ Moléculaire et La Propriété Ferromagnétique. J. Phys. Theor. Appl. 6, 661–690. doi: 10.1051/jphystap:019070060066100

CrossRef Full Text | Google Scholar

Weiss, Y. (2001). Comparing the Mean Field Method and Belief Propagation for Approximate Inference in MRFs. Burlington, MA: Morgan Kaufmann.

Google Scholar

Wens, V., Marty, B., Mary, A., Bourguignon, M., Beeck, M. O. d., et al. (2015). A geometric correction scheme for spatial leakage effects in MEG/EEG seed-based functional connectivity mapping. Hum. Brain Mapp. 36, 4604–4621. doi: 10.1002/hbm.22943

PubMed Abstract | CrossRef Full Text | Google Scholar

Wipf, D., and Nagarajan, S. (2009). A unified Bayesian framework for MEG/EEG source imaging. Neuroimage 44, 947–966. doi: 10.1016/j.neuroimage.2008.02.059

PubMed Abstract | CrossRef Full Text | Google Scholar

Wipf, D. P., Ramirez, R. R., Palmer, J. A., Makeig, S., and Rao, B. D. (2006). Automatic Relevance Determination for Source Localization with MEG and EEG Data. Technical Report. San Diego, CA: University of California.

PubMed Abstract | Google Scholar

Yedidia, J. S., Freeman, W. T., and Weiss, Y. (2003). Understanding belief propagation and its generalizations. Exploring Art. Intell. New Millennium 8, 18–9448. doi: 10.5555/779343.779352

CrossRef Full Text | Google Scholar

Yeredor, A. (2000). The joint MAP-ml criterion and its relation to ml and to extended least-squares. IEEE Trans. Signal Process. 48, 3484–3492. doi: 10.1109/78.887041

CrossRef Full Text | Google Scholar

Yuan, M., and Lin, Y. (2006). Model selection and estimation in regression with grouped variables. J. R. Stat. Soc. Ser. B 68, 49–67. doi: 10.1111/j.1467-9868.2005.00532.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Lai, Z., Xu, Y., Shao, L., Wu, J., and Xie, G.-S. (2017). Discriminative elastic-net regularized linear regression. IEEE Trans. Image Process. 26, 1466–1481. doi: 10.1109/TIP.2017.2651396

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, S., Jayasumana, S., Romera-Paredes, B., Vineet, V., Su, Z., Du, D., et al. (2015). “Conditional random fields as recurrent neural networks,” in Proceedings of the IEEE International Conference on Computer Vision (Santiago), 1529–1537. doi: 10.1109/ICCV.2015.179

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, H., and Hastie, T. (2005). Regularization and variable selection via the elastic net. J. R. Stat. Soc. Ser. B 67, 301–320. doi: 10.1111/j.1467-9868.2005.00503.x

CrossRef Full Text | Google Scholar