Abstract

In this letter, we present a new method for integration of sensor-based multifrequency bands of electroencephalography and magnetoencephalography data sets into a voxel-based structural-temporal magnetic resonance imaging analysis by utilizing the general joint estimation using entropy regularization (JESTER) framework. This allows enhancement of the spatial-temporal localization of brain function and the ability to relate it to morphological features and structural connectivity. This method has broad implications for both basic neuroscience research and clinical neuroscience focused on identifying disease-relevant biomarkers by enhancing the spatial-temporal resolution of the estimates derived from current neuroimaging modalities, thereby providing a better picture of the normal human brain in basic neuroimaging experiments and variations associated with disease states.

1  Introduction

Three cutting-edge technologies have emerged as the primary noninvasive functional neuroimaging modalities for studying the human brain: functional magnetic resonance imaging (FMRI), electroencephalography (EEG), and magnetoencephalography (MEG) (Churchland & Sejnowski, 1988). Each of these provides unique information but at different spatial and temporal scales, and thus they have their own strengths and weaknesses. The measured signal changes in FMRI are related to changes in blood oxygenation (Buxton & Frank, 1997) and are thus only indirectly related to neural activity and occur on a very slow timescale ( s). EEG and MEG capture the direct effects of neural activity with exquisite temporal resolution ( ms), but with very poor spatial resolution and can measure activity only in the cortex because electrical and magnetic fields decay rapidly with distance, whereas MRI can detect activity throughout the entire brain with excellent spatial resolution mm).

The ability to leverage the strengths of these different acquisition modalities by combining their analysis into a single framework would thus have the potential to produce estimates of structure-function relationships in the brain with greater spatial-temporal resolution than any of these methods taken individually. This potential advantage of multimodal integration of EEG, MEG, and neuroimaging with MRI (neuro-MRI) signals is well recognized and has been actively pursued in the past two decades using analysis approaches ranging from simultaneous EEG or MRI acquisition to reconstruction algorithms that model the neuronal generators of EEG and MEG potentials based on FMRI activations obtained under the same task parameters (Sommer, Meinhardt, & Volz, 2003; Singh, Barnes, Hillebrand, Forde, & Williams, 2002; Vanni, Warnking, Delon-Martin, Bullier, & Segebarth, 2004). However, current approaches typically rely on a set of restrictive assumptions, and combined with numerous conceptual and methodological challenges, the ability to fully combine modalities for the purpose of noninvasive study of human brain structure-function relationships remains an elusive capability.

This problem of combining multiple imaging modalities is ubiquitous across a wide range of scientific fields. In fact, it is an important problem in neuro-MRI, where the three primary methods—high spatial resolution anatomical data (HRA), diffusion tensor imaging (DTI), and functional MRI (either resting state FMRI (rsFMRI) or task based)—provide multivariate information on, respectively, brain morphology, neural connectivity, and brain activity at different spatial and temporal resolutions. To address this problem, we have recently developed a theoretical and computational framework for joint ESTimation using entropy regularization (JESTER) (Galinsky & Frank, 2017b) of the structural and functional parameters provided by these modalities, using each to constrain the others and thus significantly improving the veracity of the estimation procedure. The approach employs probabilistic framework to generate both intra- and inter modality coupling in a consistent manner expressed through the entropy spectrum pathways (ESP) (Frank & Galinsky, 2014). The JESTER method integrates three approaches that we developed for separately analyzing each modality:

  1. Spherical wave decomposition (SWD) (Galinsky & Frank, 2014) for shape analysis and segmentation of HRA data

  2. Geometrical optics guided by entropy spectrum pathways (GO-ESP) (Frank & Galinsky, 2014; Galinsky & Frank, 2015) for simultaneous estimation of local diffusion and global tractography from DTI data

  3. Entropy field decomposition EFD (Frank & Galinsky, 2016a, 2016b), for characterization and estimation of space-time brain activation patterns from rsFMRI, in conjunction with a new nonlinear registration method (SYM-REG) (Galinsky & Frank, 2017a, 2017b) capable of robustly and efficiently combining multimodal and multisubject data.

From a pure terminology perspective, in this letter, the EFD will refer to a procedure for extraction of a sequence of spatiotemporal patterns (activation modes) from single-modality space-time volume. And an extension of the EFD procedure that employs a combination of multiple MRI modalities (with the help of SWD and GO-ESP) for each spatial location will be referred to as the JESTER method.

It is important to emphasize that the JESTER framework is very general, may be used with data of different structure and extent, including scalar, vector, tensorial or time-series multidimensional data, and is not specific to MRI. In this letter we expanded its capabilities by incorporating both EEG and MEG with the three MRI modalities to enhance the spatial-temporal localization of brain function and relate it to morphological features and structural connectivity. This enhanced version of JESTER to incorporate cortical electrophysiological measurements is called CORT-JESTER.

2  Method

In order to be able to include EEG/MEG data sets in our joint estimation scheme, we first need to develop an approximation for the volumetric distribution of electrostatic potential inside the MRI domain. In a most general form, this approximation can be derived using Maxwell equations in a medium. Nevertheless, this is inherently an ill–posed inverse problem that requires regularization, and available multimodal MRI volumes can be used as constraints to restrict a space of available solutions.

2.1  Theory

Brain electromagnetic activity can be described using Gauss's and Ampere's laws from macroscopic Maxwell equations,
formula
where is the displacement field and is the magnetizing field. The source terms in these equations include a density of free charges and a free current density . We will assume linear electrostatic and magnetostatic properties of the brain tissues, that is, use the relations , , where is the electric field, is the magnetic field, and and are the permittivity and permeability coefficients, which may depend on the frequency .
Alternatively, both of these equations can be expressed in the form of a single charge continuity equation:
formula
Using electrostatic potential () and Ohm's law , the charge continuity equation can be rewritten as
formula
2.1
where is a local tissue conductivity tensor. Taking the temporal Fourier transform (i.e., replacing , ), this can be written in tensor form as
formula
2.2
where a summation is assumed over repeated indices.
Rearranging the above expression to separate the isotropic and homogeneous terms gives
formula
2.3
written in terms of the operators
formula
where is an isotropic local conductivity . Terms in 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 electrostatic or magnetostatic signal detected by the EEG/MEG sensors. The first term () corresponds to areas with sudden change in permittivity (e.g., the white-gray matter interface). The second term () corresponds to regions where the conductivity gradient is the strongest (the gray matter/CSF boundary). Finally, the last term () includes areas with the strongest conductivity anisotropies (e.g., input from major white matter tracts).

2.2  Inverse Solution

For this letter, it is assumed that we have the three standard neuro-MRI acquisitions (HRA, DTI, rsFMRI). The data from these modalities can be used to constrain the values of and . The DTI data allow the construction of estimates of the conductivity tensor anisotropy, whereas HRA and rsFMRI data are useful for tissue segmentation and assignment of mean values for permittivity and conductivity (see Table 1).

Table 1:
Dielectric Properties for Brain Tissues at 10 Hz Frequency.
TissueSourcePermittivityConductivity (S/m)
Gray matter Gray matter 4.07E+7 2.75E−2 
White matter White matter 2.76E+7 2.77E−2 
Cerebrospinal fluid Cerebrospinal fluid 1.09E+2 2.00E+0 
Skull cancellous Bone cancellous 1.00E+7 7.56E−2 
Skull cortical Bone cortical 5.52E+4 2.00E−2 
Skin Skin (dry) 1.14E+3 2.00E−4 
TissueSourcePermittivityConductivity (S/m)
Gray matter Gray matter 4.07E+7 2.75E−2 
White matter White matter 2.76E+7 2.77E−2 
Cerebrospinal fluid Cerebrospinal fluid 1.09E+2 2.00E+0 
Skull cancellous Bone cancellous 1.00E+7 7.56E−2 
Skull cortical Bone cortical 5.52E+4 2.00E−2 
Skin Skin (dry) 1.14E+3 2.00E−4 

Notes: The dielectric parameters are based on the Gabriel dispersion relationships (Gabriel, Lau, & Gabriel, 1996a, 1996b). The tissue permittivities are shown normalized by the permittivity of vacuum F/m.

An approximate solution for the potential across an entire MRI brain volume can be constructed iteratively as
formula
2.4
formula
2.5
where a single iteration forward solution can be found either by using the Green function or applying a Fourier-space pseudospectral approach (Gottlieb & Orszag, 1977).

We should emphasize that the use of a pseudospectral approach provides a number of important advantages over finite/boundary element (FEM/BEM) approaches most commonly performed for electrostatic modeling of brain activity (Gramfort, Papadopoulo, Olivi, & Clerc, 2010; Kybic et al., 2005; von Ellenrieder, Valdes-Hernandez, & Muravchik, 2009; Gutierrez & Nehorai, 2008; Schimpf, Ramon, & Haueisen, 2002; Ermer, Mosher, Baillet, & Leah, 2001; Mosher, Leahy, & Lewis, 1999). The pseudospectral formulation does not use surface meshes and, as a consequence, does not require limiting the location of sources by a single surface (or small number of surfaces) with a fixed number of surface-pinned static dipole sources. This is also true for target locations, as typically FEM/BEM approaches obtain the solution only at a small number of fixed position sensors. The pseudospectral approach is able to find a time-dependent spatial distribution of electrostatic potential at every space-time location of multidimensional volume as a superposition of source inputs from every voxel of the same volume. The forward solution for the electrostatic potential can be constructed through a direct summation, using static Green function and restricting the number of sources or targets by a selected set of mesh locations, thus obtaining something equivalent to an FEM/BEM solver. Alternatively, it can be obtained as a sum of inputs from all spatial and temporal Fourier modes using frequency or wave number domain for a description of both a Green function and volume distribution of sources, resulting in pseudo-spectral formulation (where the use of fast Fourier transform permits effective numerical implementation). This approach includes the distribution of both electrostatic and geometric properties of the media (conductivity, permittivity, anisotropy, inhomogeneity) at every location throughout the volume (guided by MRI acquisition and analysis); therefore, it models wave-like signal propagation inside the volume and should be able to describe and uncover significantly more complex dynamical behavior of the sources of the electrostatic activity recorded at the sensor locations. To illustrate this and show some examples of this complex activity we included in Figures 1 and 2 hyperlinks to videos of the space-time variations resulting from some rsFMRI modes used as space-time sources.

Figure 1:

Several time shots of skull electric field distributions produced by different FMRI functional modes estimated using JESTER. top left, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-1-2-b.mp4; top right, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-1-3-b.mp4; bottom left, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-2-3-b.mp4; bottom right, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-3-4-b.mp4.

Figure 1:

Several time shots of skull electric field distributions produced by different FMRI functional modes estimated using JESTER. top left, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-1-2-b.mp4; top right, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-1-3-b.mp4; bottom left, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-2-3-b.mp4; bottom right, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-3-4-b.mp4.

Figure 2:

Several time shots of brain electric field distributions produced by the same FMRI functional modes as in Figure 1 estimated using JESTER. top left, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-1-2-i.mp4; top right, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-1-3-i.mp4; bottom left, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-2-3-i.mp4; bottomright, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-3-4-i.mp4.

Figure 2:

Several time shots of brain electric field distributions produced by the same FMRI functional modes as in Figure 1 estimated using JESTER. top left, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-1-2-i.mp4; top right, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-1-3-i.mp4; bottom left, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-2-3-i.mp4; bottomright, http://abeta.ucsd.edu/Videos.git/fmri2eeg/eeg-3-4-i.mp4.

An array of EEG sensors at locations can be used to constrain the solution and find coefficients by minimizing the deviation functional,
formula
2.6
with an additional requirement that an approximation error is bounded by a chosen accuracy ,
formula
2.7
where is a potential at frequency detected by the sensor . The iterations will be stopped if the convergence cannot be achieved (i.e., . As an initial approximation , any appropriate fit to equation 2.6 can be used that satisfies , for example, a linear function of coordinates with coefficients determined from the least square fit, equation 2.6.

This approach for finding an inverse solution can be applied to an array of MEG sensors without any modifications except a difference in a form of the deviation functional, where instead of an electrostatic potential at sensor location , either magnetic flux or the projection of its gradient should be substituted for either magnetometers or gradiometers, respectively (both quantities can be derived from free vacuum form of Maxwell equations).

2.3  Wave Activation Modes

The potential , equation 2.5, is the central quantity that will be used to derive various measures. It can be calculated over arbitrary frequency ranges , and thus it is possible to calculate it over the standard alpha, beta, and delta bands (or any parts or combinations of them) typically identified in EEG. These potentials can then be converted to the time domain and used for EFD analysis (Frank & Galinsky, 2016a, 2016b), generating time-space brain activation modes ranked by the power. Alternatively, the estimated potentials can be used in the joint estimation scheme presented in Galinsky and Frank (2017a) as an additional modality in the intermodality coupling matrix , resulting in cortically electrophysiologically enhanced version of JESTER (or CORT-JESTER) that can provide significant time resolution improvement over low-frequency FMRI data.

The details of the mathematical formulation of EEG/MEG specific coupling matrix as well as a short summary of JESTER are included in appendix C.

2.4  Computational Implementation

The multisensor EEG and multimodal MRI data were used for coupling by generation of frequency-dependent inverse Green function for every voxel inside the high-resolution MRI domain using multiple-layer classification of head tissues (GM, WM, CSF, scalp and skull; see Table 1) with the spherical wave decomposition (SWD) (Galinsky & Frank, 2014)--based segmentation of the HRA T1 volumes.

The processing of both resting-state and task-based EEG data sets HRA, DTI, and rsFMRI data sets involves the following steps:

  1. Segmentation of anatomical data sets using the SWD, including skull stripping and gray/white matter extraction

  2. Registration of DTI and rsFMRI data sets with the HRA data set and extraction of functionally active modes (regions) of the brain

  3. Generation of anisotropy conductivity tissue maps from DTI data

  4. Generation of inhomogeneous permittivity tissue maps (i.e., effectively restricting possible source areas for the EEG activity) using a combination of rsFMRI activation modes and HRA/DTI gray matter, including anisotropy estimation

  5. Generation of a frequency-dependent inverse Green function for each EEG recorded frequency and solving for (see equation 2.5) that includes MRI-based regularization constraints

  6. Generation of FMRI-like volumes of signal time courses in each voxel for various EEG bands (including a low-frequency FMRI-like range, as well as the standard alpha, beta, and delta bands)

  7. Extraction of space-time modes from FMRI-like EEG volumes using the JESTER procedure

3  Results and Discussion

3.1  Data

Several independent sources of MRI, EEG, and MEG data sets were used for testing and validation of the CORT-JESTER algorithms. The variety of selected data sets allows targeting different aspects of algorithm performance. This may include, for example, assessment of mode repeatability and both intra- and intersubject similarity for different subjects, sessions, and stimuli; assessment of importance and effects of different modalities in simultaneous EEG/rsFMRI acquisitions; and comparison of modes between EEG and MEG. Some of these data sets are publicly accessible and available from the open source Human Connectome Project (Van Essen et al., 2012, 2013; Sotiropoulos et al., 2013) and include the HRA T1 volumes, the DTI volumes, as well as resting-state and task-based FMRI and MEG acquisitions (CONNECTOME-D). The HCP data sets were collected on the customized Siemens 3T Connectom scanner, which is a modified 3T Skyra system (MAGNETOM Skyra Siemens Healthcare), housed at the MGH/HST Athinoula A. Martinos Center for Biomedical Imaging (see Setsompop et al., 2013, for details of the scanner design and implementation). A 64-channel, tight-fitting brain array coil (Keil et al., 2013) was used for data acquisition. The data set contains 96 slices of 140 140 matrix (1.5 mm linear voxel size) at four levels of diffusion sensitizations (b-values b = 1k, 3000, 5000, and 10,000 s/mm) distributed over 552 total q-vectors. HCP MEG data acquisition was performed on a whole head MAGNES 3600 (4D Neuroimaging, San Diego, CA) system housed in a magnetically shielded room, located at the Saint Louis University (SLU) medical campus. The sampling rate was selected to be as high as possible (2034.51 Hz) while collecting all channels (248 magnetometer channels together with 23 reference channels). Bandwidth was set (at DC, 400 Hz) to capture physiological signals.

Several HRA, DTI and rsFMRI data sets were collected locally at the UCSD Center for Functional MRI (CFMRI). These data were acquired with a 3T GE Discovery MR750 whole body system (CFMRI-D). The anatomical T1 volumes have 168 256 256 voxel size with 1.2 0.9375 0.9375 mm resolution. A multiband DTI EPI acquisition (Setsompop et al., 2011) developed at the CFMRI employed three simultaneous slice excitations to acquire data with three diffusion sensitizations (at b-values b = 1000/2000/3000 s/mm) for 30, 45, and 65 different diffusion gradients (respectively) uniformly distributed over a unit sphere. Several baseline (b = 0) images were also recorded. The data were reconstructed offline using CFMRI's multiband reconstruction routines. The DWI data sets have 100 100 72 matrix size with 2 2 2 mm3 resolution. Whole brain BOLD resting-state data were acquired over 30 axial slices using an echo planar imaging (EPI) sequence (flip angle = , slice thickness = 4 mm, slice gap = 1 mm, FOV = 24 cm, TE = 30 ms, TR = 1.8 s, matrix size = ). Further details are available in Wong, Olafsson, Tal, and Liu (2013). All data were preprocessed using the standard preprocessing analysis pathway at the CFMRI (as described in Wong et al., 2013). Nuisance terms were removed from the resting-state BOLD time series through multiple linear regression. These nuisance regressors included (1) linear and quadratic trends, (2) six motion parameters estimated during image coregistration and their first derivatives, (3) RETROICOR (second-order Fourier series; Glover, Li, & Ress, 2000), and RVHRCOR (Chang & Glover, 2009) physiological noise terms calculated from the cardiac and respiratory signals, and (4) the mean BOLD signals calculated from white matter and CSF regions and their first respective derivatives, where these regions were defined using partial volume thresholds of 0.99 for each tissue type and morphological erosion of two voxels in each direction to minimize partial voluming with gray matter.

For EEG and same-subject MRI CORT-JESTER, testing and validation of several data sets from two unrelated studies were used.

The first study concentrated on the detection of the effects of medication on the resting-state EEG and FMRI. All the data sets for this study were collected at the Laureate Institute for Brain Research (Zotev et al., 2016) and include simultaneously acquired EEG and rsFMRI recordings at 4 ms for 500 s as well as HRA T1 volumes (LIBR-D). In addition, these simultaneous EEG and rsFMRI acquisitions for each subject were repeated three times. The clinical results for this study will be reported elsewhere. In this letter, we used these simultaneous and repeated EEG and rsFMRI acquisitions to study the repeatability of the CORT-JESTER approach.

The second set of multimodal (EEG and MRI) data sets was acquired by the Javitt Group at the Nathan Kline Institute for Psychiatric Research (Javitt, 2009; Javitt & Sweet, 2015). The data include task and resting-state EEG and FMRI recordings at 2 ms for about 30 min total and HRA T1 volumes (NKIPR-D). The task-based EEG and FMRI data sets contain several acquisitions of response for similar stimuli that can be used for test-retest purposes. The clinical results for this study will be reported elsewhere.

3.2  FMRI Modes as EEG Sources

Functional modes (Frank & Galinsky, 2016b) estimated from rsFMRI (CFMRI-D) with JESTER (Galinsky & Frank, 2017b) show the brain's functional organization and functional connectivity through a number of consistent networks at different stages of consciousness and thus represent specific patterns of synchronous activity. This synchronism can be assumed to be related to different space-time sources of electrophysiological activity. Therefore, we can use those modes to find what low-frequency response they will generate at EEG/MEG sensor positions. Several examples of generated space-time distributions of the electric fields at the skull and the cortex are shown in Figures 1 and 2 with embedded and hyperlinked videos of the space-time variations.

3.3  MEG Data

MEG data sets (CONNECTOME-D) with a relatively large number of sensors comparing to EEG (248 magnetometer and 23 reference channels versus 32 or 64 head sensors and 1 reference) allow obtaining detailed space-time wave activity modes. Several examples of these modes are shown in Figure 3. The high-resolution anatomical and diffusion-weighted data sets were used for tissue classification and assignment of electrostatic magnetostatic properties for each voxel. Selected-at-random activation modes include bilateral anterior insula, which is among the most common activation patterns in rsFMRI, as well as the area that is most commonly dysfunctional in psychiatric disorders (Goodkind et al., 2015), right caudate, and medial frontal gyrus.

Figure 3:

Demonstration of activation modes created from MEG data using JESTER and constrained by high-resolution anatomical MRI data. Activated regions are (left) bilateral anterior insula, which is among the most common activation patterns in rsFMRI, (middle) right caudate, and (right) medial frontal gyrus.

Figure 3:

Demonstration of activation modes created from MEG data using JESTER and constrained by high-resolution anatomical MRI data. Activated regions are (left) bilateral anterior insula, which is among the most common activation patterns in rsFMRI, (middle) right caudate, and (right) medial frontal gyrus.

No rsFMRI data were used for the CONNECTOME-D analysis. Generally the combination of rsFMRI spatial activation modes and and HRA/DTI gray matter estimations is used only for setting the volumetric inhomogeneous distribution of permittivity. That is, rsFMRI modes can constrain source regions only indirectly by their influence in shaping the volumetric permittivity distribution. The additional information that rsFMRI provides is expected to result in a relatively small (maybe even marginal) overall improvement due to low resolution. Nevertheless, even without the use of rsFMRI modes, the MEG activation patterns in Figure 3 still show patterns that are among the most common activation patterns observed by rsFMRI.

3.4  EEG Data

The first set of EEG data (LIBR-D) was recorded with relatively sparse spatial array of sensors (32 head sensors plus 1 reference) with temporal resolution of 4 ms. The data acquisition was completed during a resting state without task-based stimuli simultaneously with acquisition of functional MRI volumes. For each subject, three recording sessions were conducted. The CORT-JESTER procedure was used to generate a new volume for each EEG set with the same spatial and temporal resolution as the FMRI data set ( and voxels with 237 time points). Each volume was then used in an EFD procedure to generate activation modes in the same way it is done for FMRI volumes. Figure 4 shows two randomly selected modes for two subjects for all three sessions. Both subjects and modes show very good mode similarity between sessions and, at the same time, some noticeable subject differences. To quantitatively evaluate mode repeatability, we used the intraclass correlation coefficient (ICC; Sokal & Rohlf, 1995) calculated as a ratio of variance across mode groups to the total variance. For the two mode groups shown in Figure 4, the ICC is equal to 0.9985. For several other mode groups (not shown), ICC values from 0.9 to 0.99 and higher were obtained.

Figure 4:

EEG functional modes detected from volumes generated by CORT-JESTER. The modes display excellent repeatability between sessions, with the intraclass correlation coefficient (Sokal & Rohlf, 1995) equal to 0.9985, but at the same time provide noticeable subject differences. The modes were generated with the same spatiotemporal resolution as was used for FMRI acquisition ( and voxels with 237 time points).

Figure 4:

EEG functional modes detected from volumes generated by CORT-JESTER. The modes display excellent repeatability between sessions, with the intraclass correlation coefficient (Sokal & Rohlf, 1995) equal to 0.9985, but at the same time provide noticeable subject differences. The modes were generated with the same spatiotemporal resolution as was used for FMRI acquisition ( and voxels with 237 time points).

The LIBR-D set includes simultaneously acquired EEG and FMRI volumes (both task based and resting state) and potentially allows conducting a study of correspondence between EEG and FMRI. This is clearly a very important question, and a number of studies have appeared recently that focus on simultaneous acquisition techniques and provide attempts to confirm (or rebuff) the existence of correlations between EEG and FMRI data at different frequency bands or under different acquisition conditions (Meyer, van Oort, & Barth, 2013; Musso, Brinkmeyer, Mobascher, Warbrick, & Winterer, 2010; Fellner et al., 2016; Meyer, Janssen, Van Oort, Beckmann, & Barth, 2013; Chang, Liu, Chen, Liu, & Duyn, 2013; Mantini, Perrucci, Del Gratta, Romani, & Corbetta, 2007). While we agree that the question of FMRI-EEG concordance is important, the complexity of the subject, and proving or disproving the existence of the EEG-rsFMRI concordance, is clearly beyond the scope of this paper. We emphasize that the main purpose of our letter is to devise a method for the use of complementary features from EEG and FMRI modalities (Allen, Damaraju, Eichele, Wu, & Calhoun, 2018), such as difference in temporal and spatial resolution, to enhance estimation of brain activation modes.

Figure 5 shows several EEG functional modes generated for alpha band using higher spatial resolution (the same resolution as in HRA acquisition, and voxels with 237 time points).

Figure 5:

Resting-state EEG functional modes generated for alpha band using HRA spatial resolution ( and voxels).

Figure 5:

Resting-state EEG functional modes generated for alpha band using HRA spatial resolution ( and voxels).

The second set of EEG data (NKIPR-D) was recorded with a higher number of sensors (64 head sensors plus 1 reference) and and higher temporal resolution as well (2 ms). The data acquisition included both resting state and stimuli for various tasks. Figure 6 shows several EEG functional modes generated for alpha band task stimuli (the spatial resolution from HRA volume is and voxels, and the number of time points is 123). The modes show very good similarity between similar stimuli in each subject.

Figure 6:

Task-based EEG functional modes generated by CORT-JESTER for two subjects using alpha EEG band. The first three columns show three spatial projections, and the fourth column shows a temporal course (using seconds for the duration units and arbitrary normalization for the amplitude). The modes show very good similarity between similar stimuli in each subject.

Figure 6:

Task-based EEG functional modes generated by CORT-JESTER for two subjects using alpha EEG band. The first three columns show three spatial projections, and the fourth column shows a temporal course (using seconds for the duration units and arbitrary normalization for the amplitude). The modes show very good similarity between similar stimuli in each subject.

4  Conclusion

We presented an extension of the theory for the joint estimation of the structural-functional brain modes that has been initially applied to the three primary neuro-MRI modalities—high-resolution anatomical (HRA) data, diffusion tensor imaging (DTI) data, and resting-state functional MRI (rsFMRI) data—to include high temporal resolution EEG and MEG modalities.

Joint estimation refers to estimation of activation modes that includes (as a first step) the EEG time course reconstruction at each voxel location using parameters estimated from the various modalities of MRI, followed by the multimodal activation mode estimation that includes a combination of both MRI and EEG/MEG modalities.

We applied the method to multimodal MRI/EEG/MEG resting and task–based data sets showing good repeatability and similarity between subjects, with the intraclass correlation coefficient ranging from 0.9 to higher than 0.99. This method may potentially have broad implications for both basic neuroscience and clinical studies by enhancing the spatial-temporal resolution of the estimates derived from current neuroimaging modalities, thereby providing a better picture of the normal human brain in basic neuroimaging experiments and variations associated with disease states.

Appendix A:  EFD Summary

The information Hamiltonian can be written (Enßlin, Frommert, & Kitaura, 2009) as
formula
A.1
where is essentially a normalizing constant that can be ignored, is an information propagator, is an information source, and means the complex conjugate transpose. is an interaction term (Enßlin et al., 2009),
formula
A.2
where terms describe the interaction strength.

When the source term , the linear information propagator , and the nonlinear interaction terms are all known, or at least some more or less accurate approximations can be used for their description, the IFT approach provides an effective, powerful, and mathematically elegant way to find an unknown signal by using the classical solution at the minimum of Hamiltonian () or with the help of summation methods (e.g., with the help of Feynman diagrams; Feynman, 1949; Enßlin et al., 2009).

But there is a whole class of problems where those terms are unknown and too complex for deriving effective and accurate approximations. In this case, the ESP method (Frank & Galinsky, 2014), based on the principle of maximum entropy (Jaynes, 1957a, 1957b), provides a general and effective way to introduce powerful prior information using coupling between different spatiotemporal points that is available from the data. This is accomplished by constructing a so-called coupling matrix that characterizes the relation between locations and in the data:
formula
A.3
Here the are Lagrange multipliers that describe the relations and depend on some function of the space-time locations and . The eigenvalues and eigenvectors of the coupling matrix ,
formula
A.4
then formally define the transition probability from location to location of the th mode (or path, as it is often called in the random walk theory) as
formula
A.5
For each transition matrix, equation A.5, there is a unique stationary distribution associated with each mode ,
formula
A.6
that satisfies
formula
A.7
where , associated with the largest eigenvalue , corresponds to the maximum entropy stationary distribution (Burda, Duda, Luck, & Waclaw, 2009).
The EFD approach (Frank & Galinsky, 2016a, 2016b) adds those coupling matrix priors into the information Hamiltonian, equation A.1, by expanding the signal into a Fourier expansion using as the basis functions:
formula
A.8
In this ESP basis, the information Hamiltonian, equation A.1, can be written as
formula
A.9
where matrix is the diagonal matrix , composed of the eigenvalues of the coupling matrix, and is the amplitude of th mode in the expansion of the source ,
formula
A.10
and the new interaction terms are
formula
A.11
For the nonlinear interaction terms , the EFD method again takes coupling into account through factorization of in powers of the coupling matrix,
formula
A.12
where , which results in
formula
A.13
Thus, the EFD approach provides a simple expression for the classical solution for the amplitudes ,
formula
A.14
through the eigenvalues and eigenvectors of coupling matrix (which may also include some noise corrections (Frank & Galinsky, 2016a, 2016b).

Appendix B:  Coupling for Different Modalities

For the HRA data set, we define a simple intensity-weighted nearest-neighbors coupling matrix as
formula
B.1
For DWI data, the GO-ESP procedure uses the spin density function , expressed with the help of the spherical wave decomposition as
formula
B.2
formula
B.3
and generate the symmetrized scale-dependent coupling matrix as
formula
B.4
where represents the dimensionless ratio of scales of dynamic displacement to the spatial (voxel) scales , is the spherical Bessel function of order , and is the spherical harmonic, with and being the polar and azimuthal angles of the vector , and similarly for the vector , and is the DW1 signal (see Galinsky & Frank, 2015, for more details).
For rsFMRI the EFD procedure employs the frequency dependent spatial coupling matrix as
formula
B.5
formula
B.6
Here is the temporal Fourier mode of the rsFMRI data with the frequency , is the eigenvector of that corresponds to the largest eigenvalue, and can include some function of the pair correlations taken to the th power, for example, a simple mean of the pair correlations, that is equivalent to a product of signal means () for a periodic signal,
formula
B.7
or a maximum correlation of mean subtracted signal,
formula
B.8
where is the pair correlation,
formula
B.9
The EFD procedure can then be extended to incorporate information from all of the modalities simultaneously by incorporating these single modality coupling matrices , , and in equation C.1 to generate the intermodality coupling matrix .

Appendix C:  JESTER Summary

Assuming that we have different modalities with the coupling matrices that all correspond to the same unknown signal , then we can construct an intermodality coupling matrix as the product of these 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 between any two space-time locations can be written in the general (equivalent) form as
formula
C.1
where the exponents can either be some constants or functions of data collected for different modalities:
formula
C.2
and represent, respectively, the data and the coupling matrix of the modality data set represented in the ESP basis and evaluated at locations and of a common reference domain ,
formula
C.3
where denotes a diffeomorphic mapping of the th modality from the reference domain to an acquisition space .
An EEG/MEG addition to JESTER employs the frequency dependent spatial coupling matrix , similar to rsFMRI, as
formula
C.4
formula
C.5
Here is the temporal Fourier mode of the EEG/MEG data with the frequency , is the eigenvector of that corresponds to the largest eigenvalue, and can include some function of the pair correlations taken to the th power, for example, a simple mean of the pair correlations, that is equivalent to a product of signal means () for a periodic signal,
formula
C.6
or a maximum correlation of a mean subtracted signal,
formula
C.7
where is the pair correlations:
formula
C.8

Acknowledgments

We thank Alec Wong and Tom Liu at the UCSD CFMRI for providing the resting-state data and Scott Sorg at the VA San Diego Health Care System for providing the diffusion weighted imaging data. L.R.F. and V.L.G. were supported by NSF grants DBI-1143389, DBI-1147260, EF-0850369, PHY-1201238, ACI-1440412, ACI-1550405, and NIH grant R01 MH096100 and R01 AG054049-02. Data were provided in part by the Human Connectome Project, WU-Minn Consortium (principal investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. Data collection and sharing for this project was provided by the Human Connectome Project (HCP; principal investigators: Bruce Rosen, Arthur W. Toga, and Van J. Weeden). HCP funding was provided by the National Institute of Dental and Craniofacial Research, the National Institute of Mental Health, and the National Institute of Neurological Disorders and Stroke. HCP data are disseminated by the Laboratory of Neuro Imaging at the University of Southern California.

References

Allen
,
E. A.
,
Damaraju
,
E.
,
Eichele
,
T.
,
Wu
,
L.
, &
Calhoun
,
V. D.
(
2018
).
EEG signatures of dynamic functional network connectivity states
.
Brain Topogr
,
31
(
1
),
101
116
.
Burda
,
Z.
,
Duda
,
J.
,
Luck
,
J.
, &
Waclaw
,
B.
(
2009
).
Localization of the maximal entropy random walk
.
Phys. Rev. Lett.
,
102
(
16
),
160602
.
Buxton
,
R. B.
, &
Frank
,
L. R.
(
1997
).
A model for the coupling between cerebral blood flow and oxygen metabolism during neural stimulation
.
J. Cerebr. Blood Flow Metab.
,
17
(
1
),
64
72
.
Chang
,
C.
, &
Glover
,
G. H.
(
2009
).
Effects of model-based physiological noise correction on default mode network anti-correlations and correlations
.
NeuroImage
,
47
(
4
),
1448
1459
.
Chang
,
C.
,
Liu
,
Z.
,
Chen
,
M. C.
,
Liu
,
X.
, &
Duyn
,
J. H.
(
2013
).
EEG correlates of time-varying BOLD functional connectivity
.
NeuroImage
,
72
,
227
236
.
Churchland
,
P. S.
, &
Sejnowski
,
T. J.
(
1988
).
Perspectives on cognitive neuroscience
.
Science
,
242
(
4879
),
741
745
.
Enßlin
,
T. A.
,
Frommert
,
M.
, &
Kitaura
,
F. S.
(
2009
).
Information field theory for cosmological perturbation reconstruction and nonlinear signal analysis
.
Phys. Rev. D
,
80
(
10
),
105005
.
Ermer
,
J. J.
,
Mosher
,
J. C.
,
Baillet
,
S.
, &
Leah
,
R. M.
(
2001
).
Rapidly recomputable EEG forward models for realistic head shapes
.
Phys. Med. Biol.
,
46
(
4
),
1265
1281
.
Fellner
,
M. C.
,
Volberg
,
G.
,
Mullinger
,
K. J.
,
Goldhacker
,
M.
,
Wimber
,
M.
,
Greenlee
,
M. W.
, &
Hanslmayr
,
S.
(
2016
).
Spurious correlations in simultaneous EEG-fMRI driven by in-scanner movement
.
NeuroImage
,
133
,
354
366
.
Feynman
,
R. P.
(
1949
).
Space-time approach to quantum electrodynamics
.
Phys. Rev.
,
76
(
6
),
769
789
.
Frank
,
L. R.
, &
Galinsky
,
V. L.
(
2014
).
Information pathways in a disordered lattice
.
Phys. Rev. E
,
89
(
3
),
11
.
Frank
,
L. R.
, &
Galinsky
,
V. L.
(
2016a
).
Detecting spatio-temporal modes in multivariate data by entropy field decomposition
.
J. Phys. A
,
49
,
395001
.
Frank
,
L. R.
, &
Galinsky
,
V. L.
(
2016b
).
Dynamic multi-scale modes of resting state brain activity detected by entropy field decomposition
.
Neural Comput.
,
28
(
9
),
1769
1811
.
Gabriel
,
S.
,
Lau
,
R. W.
, &
Gabriel
,
C.
(
1996a
).
The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues
.
Phys. Med. Biol.
,
41
(
11
),
2271
2293
.
Gabriel
,
S.
,
Lau
,
R. W.
, &
Gabriel
,
C.
(
1996b
).
The dielectric properties of biological tissues: II. Measurements in the frequency range 10 Hz to 20 GHz
.
Phys. Med. Biol.
,
41
(
11
),
2251
2269
.
Galinsky
,
V. L.
, &
Frank
,
L. R.
(
2014
).
Automated segmentation and shape characterization of volumetric data
.
NeuroImage
,
92
,
156
168
.
Galinsky
,
V. L.
, &
Frank
,
L. R.
(
2015
).
Simultaneous multi-scale diffusion estimation and tractography guided by entropy spectrum pathways
.
IEEE Trans. Med. Imag.
,
34
(
5
),
1177
1193
.
Galinsky
,
V. L.
, &
Frank
,
L. R.
(
2017a
).
Symplectomorphic registration with phase space regularization by entropy spectrum pathways
.
Manuscript submitted for publication
.
Galinsky
,
V. L.
, &
Frank
,
L. R.
(
2017b
).
A unified theory of neuro-MRI data shows scale-free nature of connectivity modes
.
Neural Comput.
,
29
,
1441
1467
. doi:10.1162/neco_a_0095
Glover
,
G. H.
,
Li
,
T. Q.
, &
Ress
,
D.
(
2000
).
Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR
.
Magn. Reson. Med.
,
44
(
1
),
162
167
.
Goodkind
,
M.
,
Eickhoff
,
S. B.
,
Oathes
,
D. J.
,
Jiang
,
Y.
,
Chang
,
A.
,
Jones-Hagata
,
L. B.
, …
Etkin
,
A.
(
2015
).
Identification of a common neurobiological substrate for mental illness
.
JAMA Psychiatry
,
72
(
4
),
305
315
.
Gottlieb
,
D.
, &
Orszag
,
S.
(
1977
).
Numerical analysis of spectral methods
.
Philadelphia
:
Society for Industrial and Applied Mathematics
. doi:10.1137/1.9781611970425
Gramfort
,
A.
,
Papadopoulo
,
T.
,
Olivi
,
E.
, &
Clerc
,
M.
(
2010
).
OpenMEEG: Opensource software for quasistatic bioelectromagnetics
.
Biomed. Eng. Online
,
9
,
45
.
Gutierrez
,
D.
, &
Nehorai
,
A.
(
2008
).
Array response kernels for EEG and MEG in multilayer ellipsoidal geometry
.
IEEE Trans Biomed Eng.
,
55
(
3
),
1103
1111
.
Javitt
,
D. C.
(
2009
).
When doors of perception close: Bottom-up models of disrupted cognition in schizophrenia
.
Annu. Rev. Clin. Psychol.
,
5
,
249
275
.
Javitt
,
D. C.
, &
Sweet
,
R. A.
(
2015
).
Auditory dysfunction in schizophrenia: Integrating clinical and basic features
.
Nat. Rev. Neurosci.
,
16
(
9
),
535
550
.
Jaynes
,
E.
(
1957a
).
Information theory and statistical mechanics
.
Physical Review
,
106
(
4
),
620
630
.
Jaynes
,
E.
(
1957b
).
Information theory and statistical mechanics. II
.
Physical Review
,
108
(
2
),
171
.
Keil
,
B.
,
Blau
,
J. N.
,
Biber
,
S.
,
Hoecht
,
P.
,
Tountcheva
,
V.
,
Setsompop
,
K.
, …
Wald
,
L. L.
(
2013
).
A 64-channel 3t array coil for accelerated brain MRI
.
Magn. Reson. Med.
,
70
(
1
),
248
258
.
Kybic
,
J.
,
Clerc
,
M.
,
Abboud
,
T.
,
Faugeras
,
O.
,
Keriven
,
R.
, &
Papadopoulo
,
T.
(
2005
).
A common formalism for the integral formulations of the forward EEG problem
.
IEEE Trans. Med. Imaging
,
24
(
1
),
12
28
.
Mantini
,
D.
,
Perrucci
,
M. G.
, Del
Gratta
,
C.
,
Romani
,
G. L.
, &
Corbetta
,
M.
(
2007
).
Electrophysiological signatures of resting state networks in the human brain
.
Proc. Natl. Acad. Sci. U.S.A.
,
104
(
32
),
13170
13175
.
Meyer
,
M. C.
,
Janssen
,
R. J.
,
Van Oort
,
E. S.
,
Beckmann
,
C. F.
, &
Barth
,
M.
(
2013
).
The quest for EEG power band correlation with ICA derived fMRI resting state networks
.
Front. Hum. Neurosci.
,
7
,
315
.
Meyer
,
M. C.
,
van Oort
,
E. S.
, &
Barth
,
M.
(
2013
).
Electrophysiological correlation patterns of resting state networks in single subjects: A combined EEG-fMRI study
.
Brain Topogr.
,
26
(
1
),
98
109
.
Mosher
,
J. C.
,
Leahy
,
R. M.
, &
Lewis
,
P. S.
(
1999
).
EEG and MEG: Forward solutions for inverse methods
.
IEEE Trans. Biomed. Eng.
,
46
(
3
),
245
259
.
Musso
,
F.
,
Brinkmeyer
,
J.
,
Mobascher
,
A.
,
Warbrick
,
T.
, &
Winterer
,
G.
(
2010
).
Spontaneous brain activity and EEG microstates. A novel EEG/fMRI analysis approach to explore resting-state networks
.
NeuroImage
,
52
(
4
),
1149
1161
.
Schimpf
,
P. H.
,
Ramon
,
C.
, &
Haueisen
,
J.
(
2002
).
Dipole models for the EEG and MEG
.
IEEE Trans. Biomed. Eng.
,
49
(
5
),
409
418
.
Setsompop
,
K.
,
Gagoski
,
B. A.
,
Polimeni
,
J. R.
,
Witzel
,
T.
,
Wedeen
,
V. J.
, &
Wald
,
L. L.
(
2011
).
Blipped-controlled aliasing in parallel imaging for simultaneous multislice echo planar imaging with reduced g-factor penalty
.
Magn. Res. Med.
,
67
(
5
),
1210
1224
.
Setsompop
,
K.
,
Kimmlingen
,
R.
,
Eberlein
,
E.
,
Witzel
,
T.
,
Cohen-Adad
,
J.
,
Mc- Nab
,
J. A.
, …
Wald
,
L. L.
(
2013
).
Pushing the limits of in vivo diffusion MRI for the human connectome project
.
NeuroImage
,
80
,
220
233
.
Singh
,
K. D.
,
Barnes
,
G. R.
,
Hillebrand
,
A.
,
Forde
,
E. M.
, &
Williams
,
A. L.
(
2002
).
Task-related changes in cortical synchronization are spatially coincident with the hemodynamic response
.
NeuroImage
,
16
(
1
),
103
114
.
Sokal
,
R.
, &
Rohlf
,
F.
(
1995
).
Biometry: The principles and practice of statistics in biological research
.
New York
:
Freeman
.
Sommer
,
J.
,
Meinhardt
,
J.
, &
Volz
,
H. P.
(
2003
).
Combined measurement of event-related potentials EPRs and FMRI
.
Acta Neurobiol. Exp. (Wars)
,
63
(
1
),
49
53
.
Sotiropoulos
,
S. N.
,
Moeller
,
S.
,
Jbabdi
,
S.
,
Xu
,
J.
,
Andersson
,
J. L.
,
Auerbach
,
E. J.
, …
Lenglet
,
C.
(
2013
).
Effects of image reconstruction on fiber orientation mapping from multichannel diffusion MRI: Reducing the noise floor using SENSE
.
Magn. Reson. Med.
,
70
(
6
),
1682
1689
.
Van Essen
,
D. C.
,
Smith
,
S. M.
,
Barch
,
D. M.
,
Behrens
,
T. E. J.
,
Yacoub
,
E.
,
Ugurbil
,
K.
, &
for the WU-Minn HCP Consortium
. (
2013
).
The WU-Minn Human Connectome Project: An overview
.
NeuroImage
,
80
(
C
),
62
79
.
Van Essen
,
D. C.
,
Ugurbil
,
K.
,
Auerbach
,
E.
,
Barch
,
D.
,
Behrens
,
T. E.
,
Bucholz
,
R.
, …
E.
,
Y.
(
2012
).
The Human Connectome Project: A data acquisition perspective
.
NeuroImage
,
62
(
4
),
2222
2231
.
Vanni
,
S.
,
Warnking
,
J.
,
Delon-Martin
,
C.
,
Bullier
,
J.
, &
Segebarth
,
C.
(
2004
).
Sequence of pattern onset responses in the human visual areas: An FMRI constrained VEP source analysis
.
NeuroImage
,
21
(
3
),
801
817
.
von
Ellenrieder
,
N.
,
Valdes-Hernandez
,
P. A.
, &
Muravchik
,
C. H.
(
2009
).
On the EEG/MEG forward problem solution for distributed cortical sources
.
Med. Biol. Eng. Comput.
,
47
(
10
),
1083
1091
.
Wong
,
C. W.
,
Olafsson
,
V.
,
Tal
,
O.
, &
Liu
,
T. T.
(
2013
).
The amplitude of the resting-state fMRI global signal is related to EEG vigilance measures
.
NeuroImage
,
83
,
983
990
.
Zotev
,
V.
,
Yuan
,
H.
,
Misaki
,
M.
,
Phillips
,
R.
,
Young
,
K. D.
,
Feldner
,
M. T.
, &
Bodurka
,
J.
(
2016
).
Correlation between amygdala BOLD activity and frontal EEG asymmetry during real-time fMRI neurofeedback training in patients with depression
.
Neuroimage Clin.
,
11
,
224
238
.