1 Main

The human brain communicates internally through exceedingly complex spatial and temporal patterns of electrical signals. Although these signals can be measured using electrodes placed on the surface of the scalp (electroencephalography, or EEG), the ability to reconstruct the spatial and temporal patterns within the brain has been thwarted by the complexity of the inverse problem: What time- (or frequency-) dependent volumetric electrical signals throughout the brain are consistent with the signal measured on the two-dimensional surface of the scalp (1, 2)? There is a long-standing belief that it is not possible to detect and reconstruct electrical activity in sub-cortical regions deep within the brain from EEG due to inherent limitations of “volume conduction” (3). However, this is not actually a physical limitation, but rather a consequence of the incomplete nature of the standard model used to characterize the EEG signal. Despite the obviously highly dynamical nature of the electrical activity that occurs within the very inhomogeneous and anisotropic composition of brain tissue, current EEG data analysis methods are still based on the assumption that the average tissue bioelectric properties (e.g., the average permittivity and conductivity ) are sufficient to describe the electric fields E in the brain. This leads to the approximation (4) which in turn leads to the assumption that the time depedence ∂E/∂t can be ignored in the “typical” frequency range of brain signals (5). This is the ubiquitous so-called “quasi-static” approximation (6).

In reality, it is precisely the anisotropic and inhomogeneous nature of brain tissue that must be taken into account in order to develop an accurate physical model of brain electromagnetic behavior, as we have described in our recently developed universal theory of brain waves called weakly evanescent transverse cortical waves (WETCOW) (7–9). The surprising consequence of this theory is the existence of electric field waves generated as a consequence of the complex tissue boundaries (e.g., surface waves) that permeate throughout the brain and are in precisely the frequency range of observed brain electrical activity. This theory explains the broad range of observed but seemingly disparate brain spatiotemporal electrical phenomena from extracellular spiking to cortical wave loops, all of which are predicated on the time dependence of the electric fields within the complex architecture of anisotropic and inhomeneous tissue within the brain. This theory is necessary to provide a solution to the EEG inverse problem which, as shown below, produces a reconstruction of brain electrical activity with high temporal resolution and spatial resolution that is comparable (or even exceeding that of) functional MRI.

2 Methods

The essence of our method (detailed in the Supplementary Material) is to reconstruct the electric field potential throughout the entire brain using our recently developed WETCOW wave model (7–9) constrained by MRI-defined tissue properties. This method is called SPatially resolved EEG Constrained with Tissue properties by Regularized Entropy (SPECTRE). This approach solves a dynamic time-dependent construction of Maxwell’s wave equations of electromagnetism in an inhomogeneous and anisotropic medium and thus is distinctly different than standard so-called ‘source localization’ methods, as discussed in Section 6.3.

3 Validation

Validation of any neuroimaging methods is problematic because it is not possible to directly measure brain activity at every location in the brain. Nevertheless, three methods are obvious candidates for assessment of SPECTRE’s validity.

The first is comparison with functional MRI (fMRI), the current method of choice for whole brain spatial localization of brain activity. However, association of fMRI with a “standard” for EEG is problematic because it is not measuring electrical activity, but the magnetization changes in hemoglobin as blood becomes deoxygenated during brain activity. The timescale and location of these changes can be vastly different than those produced by EEG signals. Nevertheless, its capability of spatially localizing activated brain regions merits a comparison. The most direct comparison is between fMRI and EEG data collected simultaneously, which guarantees that the brain activity measured is identical in both experiments. Such “simultaneous fMRI/EEG” experiments are not particularly common as collecting EEG data within an MRI scanner during imaging is notoriously difficult, and the MR imaging procedure significantly distorts the EEG signal. However, a recent open-source study provides such data which is sufficient for our purposes.

A more direct method for validating the ability of SPECTRE to reconstruct localized electrical activity can be constructed from intra-cranial EEG (iEEG) recordings collected during epilepsy studies. Such measurements consist of specially designed EEG sensors distributed linearly along a probe that is inserted deep within a brain that has been exposed by surgical removal of a portion of the skull. By selecting only these electrodes near the brain surface from the full array of electrodes, we can synthesize an artificial surface distribution of electrodes to mimic a standard non-invasive EEG experiment (albeit with a limited coverage of the brain). We have access to such data through an ongoing study which enabled this method of validation as well.

Lastly, a comparison with current “source localization” methods would seem to be in order (10). This comparison turns out to be the most problematic as these methods all employ a very different, and quite limited, physical model for the EEG signal, and suffer from computational limitations as well. Despite attempts to make a reasonably valid comparison, it was determined that this was not possible, as described below.

3.1 Validation with simultaneous fMRI/EEG visual task

It is notoriously difficult to get high quality EEG data in simultaneous fMRI/EEG studies as the presence of the rapidly varying magnetic fields present in an fMRI acquisition distort the EEG signal. However, one recent open-source simultaneous fMRI/EEG study of a well controlled visual task (the periodic flashing checkerboard) on multiple subjects ((11), available from the Nathan Kline Institute) provides important data to address this question.

The ability of SPECTRE to faithfully reconstruct the spatial distribution similar to fMRI is shown in Fig. 1. Importantly, this comparison was performed on data from a single subject, since brain activity patterns can vary significantly between individual and averaging over multiple subjects obscures specific spatial variations important for validation. In the top rows of Fig. 1 is shown the fMRI EFD mode that automatically detects the activation in the primary visual cortex. In the middle row are shown the SPECTRE modes reconstructed using the 2mm MNI anatomical atlas, chosen because it was closest in resolution (2mm³) to the fMRI data (~ 3mm³). The very close correspondence between the spatial patterns is evident.

Figure 1

The bottom rows in Fig. 1 clearly demonstrate one of the most compelling, and perhaps surprising, aspects of SPECTRE - its ability to reconstruct activation at spatial resolution significantly higher than fMRI. This is a consequence of the SPECTRE reconstruction being based on the solution of the propagation of electromagnetic wave through specific tissue morphologies and bioelectric properties, provided by arbitrary resolution anatomical MRI data. The finer the resolution of the MRI scans, the more details can be available for the reconstruction. This is of course dependent upon the number and distribution of the EEG sensors, but certainly holds for the standard array configurations used in this paper.

Although it is an almost universally believed notion that EEG and fMRI are complementary because EEG has excellent temporal resolution but poor spatial resolution, while fMRI has poor temporal resolution but good spatial resolution, in fact SPECTRE EEG reconstructions can achieve much higher intrinsic temporal and spatial resolution. Moreover, because there are no spatial distortions in SPECTRE, this mitigates one of the aspects of fMRI that most confounds spatial localization through signal loss and non-linear geometric distortions. This is shown in Fig.2.

Figure 2

It should be noted that the ‘simple’ periodic flickering checkerboard stimulus not only activates the primary visual cortex but activates other visual and supplementary fields as well, as is evident from the activity patterns in Fig. 1. A simple stimulus does not imply a simple activation pattern. This notion was a primary motivation for our development of the EFD method for fMRI (12). The activation mode reconstructions for both the fMRI and SPECTRE data are based on the EFD which detects complex non-linear interacting spatial-temporal modes of activity (12). Thus although the task is a ‘simple’ visual stimulation, our analysis is not expected to simply detect activity in only the visual cortex, as would be produced by a more standard regression approach (11), but in a more complex set of brain networks. Indeed, multiple EFD modes are produced, though we have only shown the one incorporating the primary visual cortex. As we have argued previously (12), EFD analysis is more sensitive than simple regression techniques to the complex brain activation patterns predicted by neuroscience, and less sensitive to erroneous identification of noise or non-independent modes than the independent component analysis (ICA) (12). Indeed, one of our observations from both the fMRI and EEG data used in this study (11) is the appearance of PFC activations associated with visual stimulation, which has been suggestive of conscious visual perception (13, 14). Addressing this question is beyond the scope of the current paper.

3.2 Validation with simultaneous fMRI/EEG attention paradigm

Simultaneous EEG/fMRI were collected from subjects within a standard clinical 3T MRI scanner (see Section 6.4.1 for details). The stimuli and paradigm are described in detail in (15). Briefly, bimodal stimuli consisting of short (~ 1s) streams of simple tones (600 and 1000 Hz) alternating at 10 Hz were delivered concurrently with phase-reversing (6Hz) checkerboard patterns presented at fixation. Participants were instructed to selectively attend to either the visual or auditory aspect of the bimodal stimulus and respond when the stream of stimuli in the attended modality ends.

SPECTRE processing was performed in the alpha band. The appearance of visual stimuli elicited a reduction of ongoing alpha (7–14Hz) activity (“event-related desynchronization”, ERD) over occipital cortex, believed to occur when cortical regions are brought “on-line” for information processing (16). As in previous studies, e.g. (17), attended visual stimuli elicited increased (more negative) amplitude of the alpha ERD compared to unattended stimuli (Fig. 3A,B). In contrast, unattended, compared to attended, visual stimuli elicited a greater reduction in ongoing spectral activity within the 5–15Hz frequency range over bilateral middle frontal cortex (Fig. 3C,D). We estimated the neural sources of these attention-related modulations of oscillatory activity across the 8–12Hz frequency band which encompassed both the occipital and frontal activities (Fig.3E). Their anatomical localization was remarkably consistent across several individuals (Fig.4).

Figure 3

Figure 4

A direct comparison of the activation maps derived from both fMRI and EEG using SPECTRE for single study within two subjects (i.e. without any average over studies or subjects) is shown in Fig.5. The comparison is made by choosing specific regions of interest defined in the MNI atlas (occipital cortex and cerebellum) and correlating the activation maps derived from EFD for fMRI and SPECTRE from EEG. Comparison of the similarity of activated regions in individual subjects is generally a non-trivial problem. This is particularly true in the current case where the spatial distortions in fMRI (and lack of them in SPECTRE) make measures such as mean-squared error difficult to interpret. Therefore computation of the correlation coefficient over a predefined atlas ROI is a reasonable conservative meausure. The high correlation coefficients between the maps Fig. 5 are therefore indicative of the consistency between the fMRI and SPECTRE results in the ROI. Note that this does not imply similarity over the entire region shown. Indeed, the SPECTRE results show enhanced sensivity to activation in regions not seen in the fMRI.

Figure 5

3.3 Validation through Intra-Cranial EEG recordings

While comparison with fMRI can validate the correct detection of activated brain regions and networks, as shown in the previous section, it cannot inform the question of correct detection of electrical signals, since fMRI is based on a completely different contrast mechanism related to blood oxygenation. A direct validation of SPECTRE’s ability to faithfully reconstruct deep electromagnetic activity is, to our knowledge, only achievable with one type of data: intracranial EEG (iEEG) recordings such as those used in medically refractory epilepsy patients for seizure onset localization where the electrodes are known to be adjacent to the site of electrical activity (18, 19). We analyzed an iEEG recording of a seizure localized in the left mesial temporal region acquired at Northwell Health, NY. All implanted electrodes are shown in Fig.6 (Top row) with each yellow dot depicting one recording contact. Comparing the SPECTRE reconstruction using all of the sensor data with one using only a subset of the data comprised of only the sensors on the surface of the brain (red dots in Fig.6 (Top row)) allows the quantitative assessment of how closely the results from a set of surface electrodes correspond to those produced by intra-cranial measurements recording signal very close to the sources. The results are shown for the alpha frequency band in Fig. 6 and reveal a very close correspondance between the SPECTRE mode reconstruction.

Figure 6

4 Investigation of the ‘reward circuit

Having validated the SPECTRE method directly with simultaneous fMRI/EEG, iEEG, and an attention paradigm, we investigated the ability of SPECTRE to faithfully reconstruct the well known neural ‘reward circuit’ that is one of the most important in understanding human cognition, emotion, and behavior (20–22) and is of great clinical significance in the understanding of addiction (23, 24), mood disorders (25, 26), and a variety of other conditions (27).

We demonstrate that SPECTRE using standard EEG data can accurately map human reward pathways akin to results previously only seen via fMRI. Indeed, fMRI results have highlighted a reward system within the brain that includes midbrain dopamine producing regions (the substan-sia nigra pars compacta, the ventral tegmental area), the ventral striatum, and multiple regions within the human prefrontal cortex (28). Other research using fMRI and source localization of EEG data suggests that the anterior cingulate cortex also plays a key role in reward processing (29). In a unifying theory, it has been proposed that all the aforementioned regions work together as a neural system for the optimization of reward driven behavioral change (i.e., reinforcement leaning; (30).)

This is of particular clinical significance because addictive behaviors have long been known to be subserved by specific brain regions operating in concert as the reward circuit (31–35). The reward circuit is involved in processing rewarding stimuli of any sort and in drug addiction, substances of abuse (e.g., amphetamine) increase dopamine release in a protracted and less regulated manner as compared to typical stimuli, resulting in synaptic plasticity, and altered functioning of this circuit over time.

For our analysis we used large gambling task dataset that includes 500 participants available for download from www.osf.io/65x4v/. The details of the dataset and an extensive analysis using standard EEG analysis methodologies are presented in (36). The relevant information from this study is presented in Supplementary Materials Section 6.4.2.

For each subject trial n = 10 power modes were calculated and summed to form the single space-time SPECTRE mode (see Section 6.1.2). Fig. 7 shows three orthogonal slices of the difference in EFD power summed over all modes between conditions, averaged over all subjects. Activation in key regions of the reward circuit, including the frontal lobes, anterior cingulate gyrus, accumbens, and amygdala are clearly evident. Strong negative activation (i.e., deactivation) is evident in several structures, including the supplementary motor cortex, and the parietal operculum cortex. Activation is also apparent in the lingual gyrus and around the calcarine fissure and, as expected, in bilateral subcortical structures.

Figure 7

In Fig.8 is shown the power per brain regions as defined by the Harvard-Oxford 2mm cortical (top) and subcortical (bottom) atlases. In the cortical regions (top), strong activation is apparent in the frontal cortex (medial, orbital, operculum), cingulate gyrus, paracingulate gyrus, and unsular cortex. Activation in the accumbens is apparent from the data in the sub-cortical atlas Fig.8 (bottom). These activated regions are consistent with the known elements of the human brain reward circuitry.

Figure 8

Images of statistical significance (p < .0001) are shown in Fig. 9. It should be noted that the determination of statistical significance with SPECTRE by ‘traditional’ methods is potentially misleading as they will tend to underestimate activation significance. The estimation of the modes in SPECTRE employs EFD (12, 37) which is a probabilistic formulation that by construction incorporates space-time neighborhood connectivity so that spatially and temporally coherent patterns (“clusters”) are more probable. Traditional methods have the option for “clustering” regions of activation post-hoc into their general class of techniques called “bootstrapping” or “permutation inference”. Cluster post-detection of an activation is incommensurate with our view of the estimation process wherein the clustering in space-time is a key component indicator of high-probability regions of space-time. Spatially and temporally coherent patterns maybe of low amplitude with apparent low significance by traditional means, but those intensities are within a mode that contains very high significance in cortical regions (e.g., Fig.9) which is predicted by the WETCOW model. Thus the entire SPECTRE mode, including the somewhat diffuse lower intensity regions, is significant.

Figure 9

5 Discussion and Conclusion

The ultimate goal of functional neuroimaging is to non-invasively detect and quantify the spatial and temporal variations in brain activity. Functional MRI (fMRI) and EEG have emerged as the most ubiquitous methods because they offer two important and complementary outcomes. fMRI can non-invasively detect complex spatial and relatively low frequency temporal patterns of activity related to local BOLD changes while EEG can directly detect electrical activity but without the capability of accurate spatial localization.

Ideally, one would measure the electrical activity of the brain at high temporal resolution, as done by EEG, but with the spatial localization capabilities of fMRI. However, there has been a long-standing belief that this is not possible, due to the ‘volume-conductance’ problem (3). In this paper we have shown that this is not an actual physical limitation, as it is always presented, but an artificial constraint stemming from model for how electromagnetic wave propagate within the brain that is overly simplified by neglecting local tissue anisotropy and inhomogeneity. Our recently developed WETCOW model incorporates these tissue characteristics into the electromagnetic field equations (i.e., Maxwell’s equations), which predicts the existence of previously undiscovered waves generated precisely as a consequence of the tissue anisotropy and inhomogeity. Unlike standard electromagnetic waves characterized by a frequencies w that are proportional to their spatial frequency (or wavenumber) k, (i.e., w ~ k) the WETCOW waves have an inverse relationship w ~ 1/k. This results in waves that can permeate throughout the brain, not necessarily along neuronal pathways, that are in precisely the spatial and temporal (frequency) range of observed brain electrical activity. This model provides the physical model upon which the SPECTRE is based, using tissue models derived from standard high resolution anatomical MRI data to reconstruct the modes of brain electrical activity throughout the entire brain volume.

The SPECTRE reconstruction of EEG data provides obvious significant advantages over fMRI in temporal resolution, since EEG data has very high intrinsic temporal resolution (~ 1ms) necessary to capture rapidly varying electric field variations. Moreover, the SPECTRE algorithm can specify what frequency ranges to interrogate, providing a highly flexible analysis framework for focussed investigation of particular frequency bands of interest. On the contrary, even rapid fMRI acquisition is intrinsically limited by the temporal evolution of the contrast mechanism, the BOLD signal, which is related to blood flow and thus of quite low frequency (~ 1Hz).

While the advantages of SPECTRE over fMRI in temporal resolution are clear, what is perhaps surprising is its advantages in spatial resolution. The inverse solution that estimates the electric field potential from the EEG data is based on a physical model of wave propagation from tissues whose composition and geometry are derived from high resolution anatomical MRI data. The final resolution of the SPECTRE electric field modes is that of the anatomical data which is typically significantly higher (~ .5 – 1mm)) than the resolution of an fMRI image (~ 2mm). (There are, of course, limitations depending on the number of electrodes in the EEG system.)

But it is also important to recognize that the question of resolution in fMRI is not just a question of the prescribed image resolution of the acquisition. The BOLD physical mechanism that generates the fMRI contrast is a subtle variation in the magnetic susceptibility which causes variations in the local magnetic field, that in turn alters the local signal. fMRI acquisitions are specifically designed to accentuate this effect in order to make is observable. Unfortunately, local magnetic field variations unrelated to the BOLD mechanism, in particular strong magnetic susceptibility variations due to air/tissue boundaries such as those in the sinus cavities, cause severe non-linear image distortions that effectively alter the location and shape of the affected image volume elements (voxels). This makes even the definition of ‘resolution’ problematic, as it is essentially a spatially non-linearly varying function. Such effects are absent from EEG, which is simply a set of receiving electrodes (albeit not without its own source of artifacts) (3). The SPECTRE reconstruction uses high resolution MRI data acquired with techniques specifically designed to be insensitive to these magnetic susceptibility distortions and thus of very high spatial fidelity.

In this paper, we have successfully validated the SPECTRE method using simultaneous fMRI/EEG experiments. The results not only affirmed SPECTRE’s capability to accurately reconstruct spatial distributions of neural activity from EEG data, in alignment with the concurrently acquired fMRI data, but also revealed its efficacy in identifying robust activations across subjects that were not detectable with fMRI alone. These findings underscore SPECTRE’s potential to significantly enhance the sensitivity and scope of neuroimaging analyses. Further validation was performed using intra-cranial EEG measurements from an epilepsy study with reconstruction of data from a subset of sensors on the surface of the brain were shown to be consistent with the reconstruction from all the sensors, including those directly next to the activity source. The application of SPECTRE to high resolution EEG data during a gambling task demonstrated its ability to reconstruct a well-known and important brain circuit (20–22,24,31,32,38) that has previously only been detected using fMRI. The analysis revealed significant differences in the brain networks in the alpha range 8–12Hz, consistent with previous spatially resolved fMRI experiments but the analysis is easily carried out in any user-defined frequency ranges of interest (39–42), which will be the subject of future work. The SPECTRE methodology is applicable to any EEG study and thus holds promise for a wide range of ongoing studies of basic neuroscience of reward mechanisms and in clinical applications such as addiction.

The implications for spatially resolved EEG are important not only from a scientific perspective, but from a practical perspective as well. fMRI is a much more involved and expensive procedure, requiring highly trained research or clinical applications specialists in specially designed facilities, and subjecting the subjects to a much more claustrophobic and restricted environment, with the safety concerns always present in MRI imaging experiments. On the contrary, the portability, safety, and relative ease of EEG experiments, which can be carried out in a standard research or clinical office, makes it very attractive. The high spatial and temporal resolution capabilities provided by SPECTRE to standard EEG data offer the possibility of more detailed investigations of brain activity in a wide range of both basic research and clinical settings. This method also has important implications for the democratization of medicine worldwide where there are many populations for which advanced technologies such as fMRI are prohibitive because of cost, citing issues for large specialized equipment, and lack of highly trained personnel.

Human Subjects

All participants provided informed consent approved by the University of Victoria’s Human Research Ethics Board. The iEEG data was recorded in drug-resistant epilepsy patients undergoing invasive EEG monitoring at the North Shore University Hospital (Manhasset, NY 11030, USA) for seizure onset localization. All patients provided informed written consent according to a protocol approved by the Institutional Review Board of the Feinstein Institutes for Medical Research in accordance with the Declaration of Helsinki.

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Code Availability

Software supporting this research is subject to UCSD Technology Transfer restrictions and available with appropriate agreement after reasonable request.

Additional Declarations

The authors declare no competing interests.

Acknowledgements

LRF and VLG were supported by NSF grant ACI-1550405 and NIH grant R01 AG054049. ST was supported by NIH U24 AA021695, NIH U01 AA021692, NIH U01 DA041089, NIH R01 DA057567, OK was supported by NSERC Discovery Grant RGPIN 2016-0943. AM was supported by (NIMH) R21MH123875. SB was supported by NIDCD R01DC019979. The technology in this paper is covered under US Patents 10909414, 10789713, 11270445.

6 Supplementary Material: Data and Methods

6.1 SPECTRE details

6.1.1 Solution to the inverse EEG problem

Our method for solving the inverse EEG problem can be summarized as follows. Given a standard EEG dataset from N electrodes and a high-resolution anatomical (HRA) MRI dataset with high contrast between gray matter (GM) and white matter (WM), the solution to the inverse EEG problem can be formulated as an approximation for the volumetric distribution of electrostatic potential inside the complex inhomogeneous and anisotropic tissues and complicated morphology of the MRI domain (43). A general form of this approximation is based on Maxwell’s equations in an anisotropic and inhomogeneous medium expressed as charge continuity. This is used to derive the frequency dependent electrostatic field potential ϕ that is dependent upon the electrical properties of the tissue permittivity, permeability, and conductivity. These parameters can be estimated from the HRA data. In the Fourier (i.e., frequency) domain, the electrostatic potential satisfies the equation

which is written in tensor form where a summation is assumed over repeated indices. This can be expressed in the form in terms of the operators , a frequency dependent source term and the operator

where is a local tissue conductivity tensor, is an isotropic local conductivity. Terms in square brackets show that the parts of can be interpreted in terms of different tissue characteristics and may be important for understanding the origin of sources of the electro-/magnetostatic signal detected by the EEG sensors. The first term corresponds to areas with sudden change in permittivity, e.g. the WM/GM interface. The second term corresponds to regions where the conductivity gradient is the strongest, i.e. the GM/CSF (cerebral spinal fluid) boundary. Finally, the last term includes areas with the strongest conductivity anisotropies, e.g. input from major WM tracts. The frequency and position dependent internal sources can be used to incorporate various non-linear processes including multiple frequency effects of the efficient synchronization/desynchronization by brain waves or effects of their critical dynamics. This term is ignored in the current paper.

The inverse problem can be solved by constructing an approximate solution for the potential ϕ across an entire brain volume iteratively as and (43): where a single iteration forward solution is found using a Fourier-space pseudo-spectral approach (44). Data from MRI can be used to define the complex brain tissue morphology and constrain the tissue specific values of ∑ and ε. This procedure of inverting the WETCOW brain wave model constrained by MRI-defined tissue properties is called SPatially resolved EEG Constrained with Tissue properties by Regularized Entropy (SPECTRE).

SPECTRE is flexible in its ability to incorporate relevant prior information from MRI data. HRA data is useful for tissue segmentation and assignment of mean values for permittivity and conductivity. Diffusion MRI (dMRI) data further allows construction of estimates of the conductivity tensor anisotropy. In the present study, only HRA data were used for the estimation procedure. One important practical point is that very often HRA data are not acquired for participants in EEG studies. In this case it is sufficient to use HRA data from a standard atlas, and spatially register the EEG data to the atlas. Indeed, this was case in the present paper, where the standard 2mm, 1mm, and 0.7mm resolution T1-weighted anatomical MRI Montreal Neurological Institute (MNI) (45) atlases were used, to which the EEG data were aligned using our non-linear registration algorithm SYMREG (46).

6.1.2 Mode reconstruction

The reconstructed volumetric time series of the estimated electrostatic potential field can be thought of as similar to the EM equivalent of an fMRI dataset, containing a multitude of correlated spatiotemporal patterns or “modes” of the system.

The problem then becomes one of detecting the multiple modes in complex non-linear systems. We have addressed this problem previously in our development of the entropy field decomposition (EFD) method, which is a probabilistic framework for estimating spatial-temporal modes of complex non-linear systems containing multivariate interacting fields. It is formally based on a field-theoretic mathematical formulation of Bayes’ Theorem that enables the hierarchy of multiple orders of field interactions including coupling between fields (12, 37). Its practical utility is enabled by incorporation of the theory of entropy spectrum pathways (ESP) (47), which uses the space-time correlations in each individual dataset to automatically select the very limited number of highly relevant field interactions. In short, it selects the configurations with maximum path entropy, summarized in the equilibrium (i.e., long time) distribution μ^*.

Our terminology is as follows (see Section 6.1.3 details). The k’th EFD mode of the coupled parameters α = {α₁,…, α_m} is denoted (x,t)where k = 1,…, n for some user defined n. While each of these modes provides unique information on coherent spatio-temporal activity, for characterizing the total brain activity it is often most useful and efficient to sum these modes: (x,t). In this study we utilize both individual modes and the summed modes, as appropriate.

A strength of the EFD method is that it uses prior information contained in individual datasets - there are no training datasets or averages across datasets - just the prior information contained within the single dataset of interest. This method has shown utility in resting state fMRI data (12) and in meteorology in the application to severe local storms, in particular tornadic supercells (48). The fact that this method uses prior information embedded within single datasets without the need for any ‘training’ is of significance to clinical studies in which important individual variations can be lost in the averaging process. It is also particularly important in the current paper where our validation necessitates comparison with single subject studies.

6.1.3 Entropy Field Decomposition

The entropy field decomposition (EFD) is a general probabilistic method for the estimation of spatial-temporal modes of complex non-linear, non-period, non-Gaussian multivariate data (12, 37). The goal of EFD is to estimate the field Ψ(x,t) that describes a continuous (in both space x and time t) parameter space from which the signal s_l are discrete samples of s_l = . In general, this can be done by constructing the posterior distribution of Ψ(x,t) given the data d and any prior information I that is available, via Bayes’ Theorem re-expressed in the language of field theory (49) as

where is the partition function (considered a constant in this application) and H(s,Ψ) = - ln p(s,Ψ|I) is the information Hamiltonian which takes the form

where H₀ is essentially a normalizing constant that can be ignored, j and D are the information source and propagator, and † means the complex conjugate transpose. H_i is an interaction term (49)

where terms describe the interaction strength. In highly complex non-linear systems, j,D, and are often unknown and too complex for deriving effective and accurate approximations. In this case the ESP method (47), based on the principal of maximum entropy (50,51), provides a general and effective way to introduce powerful prior information to find the most significant contributions to H(d,θ) by using coupling between different spatio-temporal points that is available from the data itself. This is accomplished by constructing a coupling matrix that characterizes the relation between locations i and j in the data Q_ij = e^-γij where the γ_ij are Lagrange multipliers that describe the relations and depend on some function of the space-time locations i and j. The eigenvalues λ_k and eigenvectors ϕ^(k) of the coupling matrix Q

then define the transition probability from location j to location i of the k’th mode as

For each transition matrix (6) there is a unique stationary distribution associated with each path k:

where μ⁽¹⁾, associated with the largest eigenvalue λ₁, corresponds to the maximum entropy stationary distribution. Note that (7) is written to emphasize that the squaring operation is performed on a pixel-wise basis.

The essence of the EFD approach (12,37) is to incorporate these coupling matrix priors into the information Hamiltonian (3) by expanding the signal s(x,t) into a Fourier expansion using {ϕ^(k)} as the basis functions

which allows expressing the information Hamiltonian (3) in this ESP basis as

where matrix Λ is the diagonal matrix Diag{λ₁,…, λ_K}, composed of the eigenvalues of the coupling matrix, and is the amplitude of kth mode in the expansion of the source j and the new interaction terms are

from which can be derived (12, 37) a simple expression for the solution for the amplitudes a_k

through the eigenvalues and eigenvectors of coupling matrix.

The EFD methods can be extended to multiple modalities by incorporating coupling between different parameters, which we call Joint Estimation with Entropy Regularization (JESTER) (52). For m = 1, …, M different modalities d^(m) with the coupling matrices Q^(m) that all correspond to the same unknown signal s, intermodality coupling matrix can be constructed as the product of the coupling matrices for the individual modalities expressed in the ESP basis and registered to a common reference frame, which we denote : That is, the joint coupling matrix is . More specifically, the joint coupling matrix L_ij between any two space-time locations (i,j) can be written in the general (equivalent) form as

where the exponents β^(m) can either be some constants or functions of data collected for different modalities where and represent, respectively, the data and the coupling matrix of the modality dataset m represented in the ESP basis and evaluated at locations r_i and r_j of a common reference domain R:

where Ψ^m : R → X denotes a diffeomorphic mapping of m-th modality from the reference domain R to an acquisition space X.

6.2 Comparison with functional MRI

Functional MRI (fMRI) has become the de facto neuroimaging method for spatial and temporal localization of brain activity. The contrast mechanism that forms the basis of fMRI is the blood oxygenation level dependent (BOLD) variations in the magnetic state of hemoglobin and its influence on the local MRI signal as a function of the local metabolism and hemodynamics (53). Consequently, the spatial and temporal characteristics of the fMRI signal are related to blood flow and metabolic dynamics, rather than direct measures of electrical activity. In particular, the signal variations will be spatially localized in vascular pathways and the temporal variations, being related to blood flow effects, are very slow compared to electrical activity. In short, the spatial-temporal dynamics measured by fMRI need not (and, in fact, will not) correspond exactly to the spatial-temporal patterns of electrical activity. Numerous experimental realities also make fMRI problematic as a gold standard. In particular, fMRI is facilitated by enhancing the sensitivity of MRI to the BOLD contrast mechanism, which requires enhancing the sensitivity to local magnetic field variations through the use of T₂-weighted pulse sequences (54) which lead to increased geometric distortions, compromising not only spatial resolution but confounding the spatial localization of the activity in a complex, non-linear fashion. Gross distortions can lead to significantly reduced signal-to-noise and even completely unrecoverable signal loss, particularly in regions near air/tissue interfaces, such as in the prefrontal cortex (PFC). Moreover, the complex non-linear interactions between the magnetic fields and physiological variations such as respiration and cardiac pulsations produce a variety of complex spatiotemporal signal distortions (55). While mitigating these artifacts is an area of very active research, they remain a serious problem for fMRI.

Nevertheless, certain very simple task-based fMRI experimental paradigms, such as finger tapping or rapidly flickering checkerboard stimuli, repeated at periodic on/off intervals, have been established as experiments that produce repeatable robust activations in known brain networks, and are commonly used as basic testbeds for assessment of analysis algorithms. When combined with simultaneous EEG acquisition, such experiments provide two different types of data that can be compared as a form of validation, with the proviso that these two methods are imaging different physical quantities.

6.3 Comparison with state-of-the-art Source Localization methods

There is a long history of attempts to spatially localize EEG activity and these are generally called “source localization” methods (56–59). These methods are fundamentally different from the SPECTRE approach as they involve numerous stringent assumptions about brain electrical activity such as a fixed set of static dipole sources, an idealized geometric model of the head reduced to a few (typically 3) shells, that spatially close points are more likely synchronized, and the smoothness of the solution. (see (60) and references therein).

These methods all implicitly assume the “quasi-static” approximation to the EM field equation which entails ignoring the time dependent terms in Maxwell’s equations, which are dependent on tissue conductivity properties which are themselves frequency dependent. The resulting solutions are therefore static, have no frequency dependence, and are insensitive to the detailed spatially variable electrical properties of the tissues. However, as discussed in detail in the development of the WETCOW model (7–9), these assumptions are incompatible with the basic physics of brain electrical activity. The SPECTRE approach is to employ the WETCOW model and solve the actual physical problem of the complete Maxwell’s equations in an inhomogeneous and anisotropic medium. The WETCOW theory provides a comprehensive framework for characterizing the propagation of EM fields through the complex brain tissue microstructure and larger scale morphology (e.g., cortical folding), and provides the dynamic solution to the electric potential field necessary to solve the EEG inverse problem.

The pseudo-spectral computational approach used in SPECTRE has some important advantages over the finite/boundary element approaches typically used for electrostatic modeling of brain activity (61–67). It does not use surface meshes and so does not require limiting the location of sources to a small number of surfaces with fixed number of static dipole sources constrained to the surfaces. And the distribution of both electrostatic and geometric properties of the media (conductivity, permittivity, anisotropy, inhomogeneity - derived from the MRI data) are incorporated at every location throughout the volume. It is thus able to find a time dependent spatial distribution of the electrostatic potential at every space-time location of a multidimensional volume as a superposition of source inputs from every voxel of the same volume (43). These traits allow it to model wave-like signal propagation inside the volume and can detect and characterize significantly more complex dynamical behavior of the sources of the electrostatic activity recorded at the sensor locations than traditional methods.

To understand intuitively why SPECTRE is capable of reconstructing EM activity through the entire brain, including deep within subcortical structures, a simple idealized example is helpful. Consider two point current sources of different frequencies, one in the cortical layer close to the scalp, the second deep within the subcortical structures of the brain. Consider a single sensor placed on the scalp collinear with the two sources. Standard source localization methods will not see the deep source, since there is no frequency dependence, and the signal falloff is simply a function of the distance from the sensor. Therefore, the close source completely dominates the signal model. Since all tomographic imaging methods (e.g., MRI, CT, etc.) depend strongly on both the spatial and temporal sampling of the measured physical system, this effective invisibility of currents in the standard quasi-static model essentially precludes the solution of the true inverse EEG problem and necessitates the artificial construction of assumed dipole distribution on pre-chosen artificial internal structures. In contrast, in SPECTRE the sources are not dipoles, but frequency sources that extend throughout the entire brain volume subject to the boundary conditions imposed by both the tissues geometry and its spatially and frequency dependent properties. The surface electrodes are assumed to be sensing EM waves emanating from the entire brain across a broad frequency spectrum limited only by the sensors. Used in conjunction with an HRA MRI data that provides the spatial distribution of the frequencydependent tissue electrical properties that constrain the possible solution, SPECTRE can invert the wave equations to provide an estimate of the spatiotemporal distribution of the electric field potential.

The problem of spatially localizating the EEG signal involves estimating the most probable distribution of electric field amplitudes given an array of sensors. This is essentially a problem of correctly modeling the physics of how electromagnetic waves propagate through the complex environment of a the convoluted brain tissue morphology and the anisotropic and inhomogeneous nature of brain tissue. The current state-of-the-art approach to this problem, called ‘source localization’ such as low resolution electromagnetic tomography or LORETA algorithm with its many variations (56–59), also called ‘EEG source imaging’ (68,69)) involves using a pre-defined brain atlas, arbitrarily placing dipole source on the surface, and calculating the contribution from these sources. Some methods propose using fMRI as a prior, which has the disadvantage of requiring fMRI acquisitions (70–73). The current source localization methods are based on a static model for the electric field caused by a fixed set of pre-defined dipole sources. (see (2) for a review of current methods). This model is inherently limited because in reality the brain’s electrical field variations are time dependent and generated by an essentially continuous distribution of sources throught the entire brain. This description is the essence of the WETCOW theory (7–9) which describes how highly coherent localized electric field phenomena, such as cortical wave loops and synchronized spiking, are produced by the complex non-linear interactions of waves across multiple spatial and temporal scales. In a typical application of the SPECTRE method, we use MNI volumetric grid with 2mm (902629 voxels), 1mm (7221032 voxels), or 0.73mm (11393280 voxels). All voxels in our models are considered sources of electromagnetic activity consistent with the local intravoxel tissue characteristics (via (1)) rather than an assumed dipolar form. In an attempt to compare SPECTRE with current state-of-the-art source localization method, we downloaded the currently available LoretaKey1 software (59), which uses a set of 6239 fixed dipoles. In order to make a fair comparison, we tried to use LoretaKey1 with a number of dipoles comparable to our 2mm number of voxels. We found that LoretaKey1 is not able to handle this size due to “out-of-memory” crashes. We scaled down the number of dipoles in 2, 4 , 10 and 20 times. Only with around 45K dipoles were we able to make Loreta run. It ran for approximately 24 hours, but then again crashed due to out-of-memory problems. Our processing with 2mm, requires around 650Mb of memory. At this point it was decided that it was not possible for Loreta to provide a result that would usefully inform the efficacy of the SPECTRE method. Our highest 0.73 resolution processing can be completed on a modern workstation taking 16–20Gb of memory in a matter of minutes.

6.4 Data

6.4.1 Attention paradigm data

Simultaneous EEG/fMRI acquisition:

Functional and structural MRI images were acquired on a Siemens 3T TIM-Trio scanner (NKI Center for Biomedical Imaging and Neuromodulation) equipped with a 32-channel phased array head coil. Structural T1 and T2 scans were collected using standard sequences. Whole-brain BOLD data was acquired with a gradient-echo EPI sequence (TR=2000ms;TE=30ms;flip angle=80°). EEG data were acquired concurrently with fMRI using an MR compatible EEG amplifier (BrainVision MR series, Brain Products, Munich, Germany) and a 64-channel MR-compatible ring electrode cap with 10–20 International System electrode placement cap. EEG data was sampled at a rate of 5 KHz EEG data were acquired at a rate of 5 KHz using BrainVision Recorder software (Brain Products). Electrocardiographic data were captured from electrodes on the backs of subjects. The reference electrode was positioned between Fz and Cz. Scanner and heartbeat artifacts were removed offline from the EEG signal using an average template subtraction procedure (74) and the data was resampled to 250Hz.

Traditional EEG analysis:

The single-trial EEG signal from each electrode was convolved with a 3-cycle Morlet wavelet computed over a 3 second window centered at the onset of each stimulus and averaged separately for each stimulus type. The averaged spectral amplitude at each time point was then baseline corrected by subtracting the mean spectral amplitude over the -200 to -50 pre-stimulus interval. Further details of post-processing and time-frequency analyses methods are described in (75–77).

6.4.2 Reward circuit data

Task and Design Acquisition

Participants completed a simple gambling task (36). On each trial, they saw a black fixation cross for 500 ms, followed by two colored squares for 500 ms, and then, the fixation cross turned gray (go cue) and participants were to select one of the two squares (square locations—left, right—were randomized on each trial) within a 2,000 ms time limit. They were then presented with a black fixation cross for 300 to 500 ms, and then, simple feedback as to their performance (“WIN” for gain, “LOSE” for loss) for 1,000 ms in black font. If the participants responded before the go cue they were instead delivered “TOO FAST” feedback and if they did not respond before the 2,000 ms time limit, it would be considered a loss. The goal of the participants was to accumulate wins by determining which of the two squares would more often lead to gains (60% vs. 10%). In this task, participants accumulated wins; however, were not paid money. They would see the same pair of colors for one block of 20 trials. They conducted six blocks of unique color pairs.

Participants

Five hundred undergraduate students were included and were recruited via the University of Victoria psychology participant pool (see (36) for details). The data was collected until 500 participants became available that were not characterized by one of the following a priori criteria: trial count after artifact rejection were less than 15 per condition, total artifact rejection exceeded 40% of trials rejected, FCz (electrode of interest) specific artifact rejection exceeded 40% of trials rejected, or independent component analysis based blink correction failed. These criteria were extremely strict to ensure clean data in the analyses, and as such a total of 637 participants were analyzed before reaching the goal of 500 clean participants. All participants had normal or corrected-to-normal vision and volunteered to take part in the experiment for extra course credit in a psychology course. All participants provided informed consent approved by the University of Victoria’s Human Research Ethics Board.

Data Acquisition and Preprocessing

EEG data were recorded from either a 64 or 32 electrode (Ag/AgCl) EEG system (ActiCAP, Brain Products, GmbH, Munich, Germany) using Brain Vision Recorder (Version 1.21.0004, Brain Products GmbH, Munich, Germany). During recording, electrodes were referenced to a common ground, impedances were kept below 20 kQ on average, data were sampled at 500 Hz, and an antialiasing low-pass filter of 245 Hz was applied via an ActiCHamp amplifier (Revision 2, Brain Products GmbH, Munich, Germany). Stimuli and EEG markers were temporally synced using a DataPixx stimulus synchronization unit (VPixx, Vision Science Solutions, Quebec, Canada)

Data were re-referenced to an average mastoid reference and filtered using a 0.1 to 30 Hz passband (Butterworth, order 4) and a 60 Hz notch filter. Correction for eye blinks was performed using EEGLAB’s independent component analysis (ICA). Components reflective of blinks were manually identified and removed via topographic maps and component loadings, and data were reconstructed. Data were then segmented from -500 to 1,500ms relative to feedback stimulus onset, baseline corrected using a -200 to 0ms window, and run through artifact rejection with 10μV/ms gradient and 100/μV maximum-minimum criteria. Data were pre-processed to identify noisy or damaged electrodes using artifact rejection trial removal rates for each electrode.

The 1 second of recorded sequence for each “WIN” or “LOSE” event were extracted from recordings for each participants (with 22 ms of pre-event sample and 488 ms of post-event sample) and combined together to form separate winning and loosing datasets. Each of those datasets were processed using SPECTRE to construct the approximate inverse solution for the potential ϕ across an entire 2 mm MNI brain volume.