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 ϱ ∈ ℝ3 (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).
Where Lvι 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) Lvι 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 . An inverse operator projects the data v(t) to explain approximately its source and produce its estimator .
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 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 . Considering that solving an inverse problem in the frequency domain should be (ideally) optimizing the inverse operator Eq. (3) from the data Fourier transform v(f) (and not from the process v(t)).
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 , 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 . In practice, this later quantity must be substituted by its estimator Eq. (4).
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 (vm (t); ∀ m). Thus, we worked with instances (vm (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).
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 () are the spectra of source activity, the Cortical Spectral Topography (CST). Off-diagonal elements are the cross-spectra reflecting functional connectivity. Optimal inverse operators 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 of the source cross-spectrum Σιι(f), may be obtained from three different types of primary information: v(t), or v(f), or . Initially, it might seem that an inverse solution 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 , 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 Lvι (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 that maximize the posterior probability (Eq. 6).
Where is the likelihood of the data conditional on the latent variable . In our context, can corresponding to be v(t), or v(f), or , and 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 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 Tvι, which holds the same linearity attribute as the forward operator Lvι, 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 Tvι 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 (Eq. 7). Here, we considered to be calculated from the set of random instances (ιm (f); ∀ m) of the latent vector process ι(f).
Once was assumed, the ABCs were simplified dramatically since a MAP for the matrix (Eq. 7) turns into a joint-MAP (Section Validation rationale) for the random vectors (ιm (f); ∀ m) that are implicit in (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 or (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 calculated from an inverse solution or . Minimizing CST distortions (for the estimator of the spectrum ) 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).
• 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(τ).
Here, the auto-covariance 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 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.”
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 the sampled estimator for Fourier transform instances (xm (f); ∀ m) with sample size M. Here, these instances (xm (f); ∀ m) represent the discrete Fourier transform applied to sample realizations (xm (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 .
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 Lvι (Eq. 12). Preserving the Gaussianity of the Fourier transform ι(f) is only possible under the linear forward operator Lvι and subsequent linear or quasilinear inverse operator Tιv (Marzetti et al., 2008; He et al., 2019; Nolte et al., 2020).
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 , and Sensor Spectral Topography (SST), a map of the sensor spectrum (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. 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 and of which we represent the topographic projections: (a) Cortical Spectral Topography (CST) and (b) Sensor Spectral Topography (SST) .
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 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 , and the forward operator Lvι.
Figure 2. Inverse operators from (a1) the MEG/EEG/ECoG sensor signal vt(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 ι3(t) and then compute (b2) their Fourier transformι3(f). Route 2 (a1→a2→b2→b3) first computes (a2) the data Fourier transform vm(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 v3(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 , and likewise, represent all types of latent variables in Eq. (13). A maximum a posteriori (MAP) Bayesian inverse operator (Eq. 14) achieves the maximum of a posteriori probability . We define the a posteriori as the conditional probability determined from the likelihood and a priori . Depending on an inverse operator can be non-linear or linear, computed numerically or analytically, also intractable or tractable.
The essential role of the a priori is to overcome the ill condition of the likelihood , in other words, to provide a unique solution to the inverse problem for the source variable . Here we formulate this a priori as a Gibbs probability density (Eq. 15).
Where α is the Gibbs temperature parameter. The Gibbs energy function 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 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 .
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 under the −log transformation. However, the Bayesian formalism is more general than the Tikhonov regularization since the Gibbs energy (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.
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 . 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 .
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).
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 Lvιι(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 Lvιι(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-normp) upon real-valued ι(t) or complex-valued ι(f) (Riesz, 1910; Rudin, 1970; Dunford and Schwartz, 1988; Bourbaki, 2013).
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 H2, 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).
Representing the most common quasilinear cases the following inverse operator , in brief notation (Eq. 20), the constrained generalized inverse of the forward operator Lιv that incorporates the regularization matrices A(f) and B(f).
Where arises from the route 2 MAP1 (Eq. 17) that incorporates the Gibbs energy H2, 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 and covariance matrix , that follows from the conjugated relation between the likelihood and a priori probabilities in the route 2 MAP1.
The cESI estimator is then (Eq. 21) for any set of inverse solution instances , or more compactly from an estimate of the data cross-spectrum (Eq. 22).
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) = Aac 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) , . 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 based on the second-type MAP (MAP2) (Eq. 23).
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 leads to a posteriori . Such a likelihood is the Complex Wishart probability density Wℂ (Eq. 11) now defined for (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).
Depending on the likelihood and a priori probability (Eq. 24), the optimal inverse operator (Eq. 23) is often non-linear and numerically intractable. In this case, which we followed in this article, optimizing requires Approximated Bayesian Computation (ABC) (Csilléry et al., 2010). The ABC we describe here employed , a non-convex but numerically tractable successive approximation to . In turn, the tractable followed from , the non-convex relaxation of the a posteriori 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 , 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 (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 and source variable by the instances and .
In our case, we considered the above p(vm(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).
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 .
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 (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 , not the more general case of Gibbs energy defined upon the instances H(ιm (f); ∀ m).
We leverage H1 (Eq. 28) the structured p = 1, q = 2-norm (square) that performs sparse regularization of the cross-spectrum topographic projection or spectrum . H1 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 .
In addition, we introduce H2 (Eq. 29) the structured p = 2, q = 2-norm (square), Shatten p = 1-norm or nuclear norm , to compensate sparse bias of H1 (Eq. 28) and regularize eigenvalues for the cross-spectrum . H2 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).
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 H1 (Eq. 28) and the nuclear norm H2 (Eq. 29) leads to the following joint a priori probability p(ιm (f); ∀ m) (Eq. 30).
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 with the square Schatten p = 2-norm (Frobenius norm) (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 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 . Solving the gamma-MAP under cross-spectral Gibbs energy 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 (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 (Eq. 31) as represented by separable (in most cases also indistinguishable) energy terms Hg(ιm (g , f); ∀ m) for each source. A self-consistent energy form Hg(ιm (g , f); ∀ m) also summarizes the interaction field g ↔ ∀g′, as a property of the original Gibbs energy H(ιm (f); ∀ m).
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 (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 , as defined by the norms considered here (Elastic nuclear quasinorm (Eq. 30) or others), fulfills such a property.
Then, if the joint a priori probability p(ιm (f); ∀ m) (Eq. 27) is written in terms of the separable Gibbs energy (Eq. 32) it turns out the Markov random field property (Eq. 33) that is summarized by 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.
Hence, we target VB—the factorizable (separable) approximation (Eq. 34) of the joint a priori p(ιm (f); ∀ m) (Eq. 33). Essential to obtain 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).
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 (Eq. 35). That is to minimize the Kullback-Leibler divergence of the approximation hqregarding to the joint a priori p (Eq. 34).
Minimizing the , 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 DKL 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 (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).
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 (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 , or directly the cross-spectrum , due to the quasilinear inverse operator (37). Quasilinear was then equivalent to iterated MAP1s for the Fourier transform (Eq. 20), positing complex-valued Gaussian likelihood p(vm(f)|ιm(f)) (Eq. 26) and a priori h(ιm(f)|γ(f)) (36).
These iterated MAP1s are for an equivalent a posteriori (Eq. 38), specified by the mean and the covariance , 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.
Targeting the variational vector γ(f) was under the marginal or hyper-parametrized likelihood for data p(vm(f)|γ(f)) (Eq. 39). This likelihood was due to integration (expectation) of p(vm(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).
However, the actual marginal likelihood p(vm(f)|γ(f)) (Eq. 39) was intractable, with approximations via the expectation of the log joint probability p(vm(f), ιm(f)|γ(f)) (Dempster et al., 1977; Liu and Rubin, 1994) under the iterated a posteriori (Eq. 40). Then an approximated marginal likelihood depended on the variational hyper-parameters (spectrum) γ(k)(f)in the current iteration.
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 was the tractable successive approximations, with the iterated a posteriori and likelihood .
Here we implemented the Bayesian learning algorithm that effectuates the gamma-MAP γ(k + 1)(f) (41) and computes the quasilinear (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 .
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. 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) for EEG sensors (green dots) and (a2) ) for MEG (ECoG) sensors (blue dots). Using the corresponding cross-spectrum data and Lead Fields (a3) and (), 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 , based on the tested method “M,” and MEG (ECoG) , based on the reference method sMNE. Employing these inverse operators, we obtain estimators for (b1) the EEG tested CST given by method “M” and (b2) σ0 the MEG (ECoG) sMNE reference CST. Incongruence between (b1) the EEG tested CST and (b2) the MEG (ECoG) reference CST is measured through (b3) eM the Earth Movers’ Distance (EMD) and cM the correlation distance (1-CORR). Using (a1) the EEG cross-spectrum data and (b3) Lead Field for EEG, we obtain the resolution operator RM. 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) rM, the Generalized Point Spread Function (GPSF) and bM 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. 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) the Earth Movers’ Distance (EMD) and 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) rM, the Generalized Point Spread Function (GPSF) and bM the Blurring for the GPSF. Elements (a1, a2) of this figure, are freely available in http://www.www.neurotycho.org/.
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 (σ0) 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 σ0 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 RM = Tιv, MLvι. Here Tιv, M is the inverse operator for the method M, and Lvι the lead field. Consider a point source uo indexed by g0, any column RM(:, g0) is then the voltages v0 produced by this source.
• The Generalized Point Spread Function (Grova et al., 2006b; Haufe et al., 2008) (GPSF) rM 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 σ0 and are the set G0. The GPSF, for any point g ∈ G in the set of cortical points, is then:
Where |𝔾0| denotes the number of elements in that set.
• The blurring measure for images (BLUR) bM 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 bM = 0 in Figure 3c1, there would be an exact coincidence between the GPSF rM non-zero values in the colormap and the blue dots in Figure 3c2. Formally, bM is defined as:
Where is the spatial dispersion around the reference point g0 and
Where 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) eM 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) cM which measures deformations from the expected collinearity between pairs of spatially distributed CSTs.
These incongruence measures (EMD) eM correlation distance (1-CORR) cM combine sensitivity to leakage and localization error.
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 rsELOR and rsLCMV 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 rssSBL was correctly centered around reference points and did not extend significantly beyond these. In other words, as can be appreciated qualitatively, rssSBL minimized leakage compared to sELORETA and LCMV.
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 rM, 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 rsELOR and rsLCMV (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 rsELOR and rsLCMV. 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 σ0 being highly frequency dependent, the GPSF pattern was rseLORETA and rsLCMV 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 rssSBL 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 bM (BLUR) (Figure 6, left-column) confirms our qualitative impressions based on the GPSF rM. The radar graphs showed a decrease in leakage of ssSBL compared to sELORETA and sLCMV (bsELOR > bssSBL, bsLCMV > bssSBL). This improvement was valid for all spectral bands (delta, theta, alpha, gamma 1, and gamma 2).
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 BM, the Blurring (BLUR), and Incongruence, employing eM the Earth Movers’ Distance (EMD), and CM 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 eM are intuitively the amount of work to deform the EEG-based CST estimator to the reference CST estimator (σM → σ0). They were one to two orders larger for sELORETA and sLCMV than SSBL (esELOR >> essSBL, and esLCMV >> essSBL).
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, csELOR was in the range of around 0.7, for all frequency bands. This behavior was congruent with that observed in the colormaps for GPSF (rsELOR and rsLCMV) (Figure 5), where the maximum values appear in opposite areas to the blue dots (reference estimation).
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, σ0) 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 rssSBL and measured in bssSBL (Figure 6, right-column), the performance of ssSBL in reducing leakage was superior compared to sELORETA (bsELOR > bssSBL) and sLCMV (bsLCMV > bssSBL).
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 rM, 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 eM (Figure 6, right-column) were consistent with those of the previous experiment (left-column), confirming the improvement of ssSBL regarding sELORETA (esELOR > essSBL by a considerably narrow margin) and LCMV (esLCMV >> essSBL by a wider margin of three orders of magnitude) for all spectral bands.
The correlation distance cM 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 (csLCMV > cLORETA>cssSBL) 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. 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).
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 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.
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.
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.
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.
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.
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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2023.978527/full#supplementary-material
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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.
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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).
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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.
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