Model driven EEG/fMRI fusion of brain oscillations

Abstract

This article reviews progress and challenges in model driven EEG/fMRI fusion with a focus on brain oscillations. Fusion is the combination of both imaging modalities based on a cascade of forward models from ensemble of post-synaptic potentials (ePSP) to net primary current densities (nPCD) to EEG; and from ePSP to vasomotor feed forward signal (VFFSS) to BOLD. In absence of a model, data driven fusion creates maps of correlations between EEG and BOLD or between estimates of nPCD and VFFS. A consistent finding has been that of positive correlations between EEG alpha power and BOLD in both frontal cortices and thalamus and of negative ones for the occipital region. For model driven fusion we formulate a neural mass EEG/fMRI model coupled to a metabolic hemodynamic model. For exploratory simulations we show that the Local Linearization (LL) method for integrating stochastic differential equations is appropriate for highly nonlinear dynamics. It has been successfully applied to small and medium sized networks, reproducing the described EEG/BOLD correlations. A new LL-algebraic method allows simulations with hundreds of thousands of neural populations, with connectivities and conduction delays estimated from diffusion weighted MRI. For parameter and state estimation, Kalman filtering combined with the LL method estimates the innovations or prediction errors. From these the likelihood of models given data are obtained. The LL-innovation estimation method has been already applied to small and medium scale models. With improved Bayesian computations the practical estimation of very large scale EEG/fMRI models shall soon be possible.

Figure 1

Strategies for EEG/fMRI data analysis. A: Simplified underlying forward models (FMs) for fusion. In a given voxel neural activity generates an ensemble of postsynaptic potentials (ePSP). Along the left branch of the diagram the temporally and spatially synchronized summated PSPs of neurons with open fields produce the primary current density (PCD). This is the PCD FM. The volume conductor properties of the head transform the PCD into EEG/MEG which is the FM for this type of signal. Along the right branch of this diagram, the ePSP generates a vasomotor feed forward signal (VFFS) via its own FM, which in turn is transformed via the hemodynamic FM into the observed BOLD signal. Note that any one of the model constructs enumerated here (ePSP, PCD, VFFS, EEG/MEG, BOLD) is time dependent and, according to the type of modeling, can be defined as either a single variable or a vector of variables. B: Fusion by measuring covariation of the EEG and BOLD. In this data driven approach to EEG/fMRI fusion the EEG is considered to have the same time evolutin as the PCD which is considered as a driver for the BOLD signal. The time varying power in an EEG band is convolved with a hemodynamic response function, h(t), and then correlated with the BOLD signal (correlation denoted by a thick arrow). Since the temporal dynamics of the EEG are taken as a surrogate for the VFFS this is an asymmetrical type of fusion. C: Fusion by measuring co‐variation of the PCD and VFFS. This is also a data driven approach in which the FMs for the EEG/MEG and BOLD are inverted, by solving respectively a spatial and temporal inverse problem to yield estimates of the PCD and VFFS. These estimates are then correlated (thick arrow) to accomplish a fusion that is symmetrical in that both modalities are given equal a priori weight. D: Model Driven Fusion by estimating the ePSP from EEG and BOLD. This is a model driven approach in which simultaneous Bayesian inversion is carried out with all FMs. In practice this involves repeated simulations with tentative values of ePSPs and other parameters and then modifying them to maximize an statistical measure of fit. One possible method for this estimation is shown in Figure 4.

Figure 2

(Color) A: Fusion by correlation of EEG time varying spectra over all channels and BOLD. Data driven asymmetrical fusion as in Figure 1b where EEG spectra convolved with hemodynamic response function is correlated with BOLD signal using multi‐linear partial least squares which detects common atoms in both modalities. Top: Spatial signatures of PCD for the EEG alpha atom estimated by f inverse solution. Bottom: spatial signature of BOLD alpha atom. B: Fusion by correlation of log PCD alpha power and log VFFS. Symmetrical data driven EEEG/fMRI fusion as in Figure 1c obtained by estimating the time course of EEG source power in the alpha band using the VARETA inverse solution, estimating the VFFS by smooth deconvolution of the BOLD signal with a standard hemodynamic response, and correlating the log of these quantities at each voxel. C: Correlation of source EEG and BOLD in neural mass network of moderate size. Correlation of estimated PCDs and VFFS obtained from a medium sized simulation of realistically connected neural masses. Connectivity measures for neural masses were estimated by means of DTI images. Correlation values threshold using the local false discovery rate (fdr). Red corresponds to positive correlations, blue to negative correlations. Simulations produced by the Local Linearization (LL) integration scheme for stochastic differential equations. D: Correlation of source EEG and BOLD in large scale neural mass network (surface of cortex and thalamus. Correlation of estimated PCDs and VFFS obtained from a very large sized simulation of realistically connected neural masses. Correlations shown for the cortical surface (above) and the thalamus (below not shown to scale) with the rostral part shown to the left. The cortical surface comprised 8203 Jansen modules were placed) and the left portion of the thalamus 438 TC/RT modules. The simulation produced (not shown) similar distributions for the right brain. All modules interconnected using the connectivity matrix shown in Figure 5. Note positive (red) PCD/BOLD correlation for caudal thalamus and negative (blue) for part of the occipital cortex marked by an arrow. Simulations produced by the approximae Local Linearization (aLL) integration scheme for stochastic differential equations.

Figure 3

Neural Mass model for EEG/fMRI fusion. This model is obtained by linking the parameters of a interconnected set of N _m neural masses to the forward models for EEG and fMRI outlined in figure 1A where averaged population values of postsynaptic potentials (PSP) serve as state variables. A: State Space Model (SSM) for a Neural Mass. Each Neural Mass (numbered as i) is shown on the right part of the figure and depicted as a component which receives two external inputs: the net input PSP z _i(t) (continuous arrow), and Gaussian white noise (dashed arrow). The neural mass produces as an output the ePSPs x ${}_{\text{2}\mathit{\text{i}}\text{−1}}^{\text{ePSP}}$ which are also the state variables for this SSM. The ePSPs will feed into other neural masses. Note that z _i(t) is the sum of the ePSPs that are emitted from other connected neural masses, amplified by the synaptic contact coefficients z _i(t) = ∑${}_{\mathit{\text{j}}\text{=1}}^{\text{N}{}_{\mathit{\text{m}}}}$ c _i,j x _2j−1(t), where c _i,j indicates connections from population j to i. By convention these synaptic connectivities will be negative for when j is an inhibitory cell. On the left is shown in more detail the sequence of operations which take input to output: 1 − z _i(t) is converted into an average pulse density of action potentials p _i(t) = S(z _i(t)) by a static nonlinear sigmoid function S(v). 2‐The noise input, with mean μ_i and standard deviation σ_i, is added to the pulse density. 3 ‐ is converted into the output x ${}_{\text{2}\mathit{\text{i}}}^{\text{ePSP}}$(t) by linear convolutions with the neural PSP impulse responses h _i(t). The complete neural mass model creates two additional signals (not shown) that will act as the two components of the VFSS signal originating the BOLD signal. These are the sum of all EPSP u _E(t) = ∑${}_{\mathit{\text{i}}\text{=1,}\mathit{\text{N}}{}_{\mathit{\text{m}}}\text{;j excitatory}}$ c _i,j x2_j−1(t) and all IPSP u _I(t) = −∑${}_{\mathit{\text{i}}\text{=1,}\mathit{\text{N}}{}_{\mathit{\text{m}}}\text{;}\mathit{\text{j}}\text{inhibitory}}$ c _i,j x2_j−1(t). B: Model for a Jansen‐Rit Cortical Module. Here N _m = 3. Pyramidal cell (Pyr), stellate cell (St), and inhibitory interneuron (Inh) populations generate ePSPs, denoted respectively by {x ${}_{\mathit{\text{1}}}^{\text{ePSP}}$ (t), x ${}_{\mathit{\text{3}}}^{\text{ePSP}}$(t), x ${}_{\mathit{\text{5}}}^{\text{ePSP}}$(t)}. The output of Pyr cells, x ${}_{\mathit{\text{1}}}^{\text{ePSP}}$(t), multiplied by c _1,2, c _3,2 drives the St and Inh populations. Pyr receives feedback from St and Inh with coefficients c _2,1, c _2,3 respectively. The only dynamic noise input is to the St Population. On the one hand, the nPCD is the transmembrane PSP of the Pyr population which is equal to the net input PSP z ₂(t) generated by the St and Inh populations z ₂(t) = c _2,1 x ${}_{\text{1}}^{\text{ePSP}}$(t) + c _2,3 x ${}_{\text{3}}^{\text{ePSP}}$(t) (EPSP‐IPSP). This is also equal to the EEG with the trivial lead field K = 1. Loosely speaking, this is as though we were actually measuring a Local Field Potential (LFP). On the other hand, The VFFS has two components , the sum of excitatory PSP and u _I(t) = c _2,3 x ${}_{\text{5}}^{\text{ePSP}}$(t) the inhibitory PSP. The VFSS components are fed into the Metabolic Hemodynamic Model (MHM) which we shall not detail here (see Appendix C) a system of ODE which depends on additional state variables x ₇(t),…,x ₁₄(t). The observed BOLD is generated by the Balloon model then is transformed to BOLD by the equation y ₂(t) = g ^BOLD (x ₁₃(t),x ₁₄(t)).

Figure 4

Local linearization (LL)‐innovation approach to estimating states and parameters of nonlinear random systems. A state space model of a dynamical system is the combination of: (a) A continuous time stochastic differential equation describing the evolution of system states x, and (b) A discrete time observation equation that explains how the observations y _t are obtained from the states. Discretization of the state equation by means of the Local Linearization scheme allows the application of sequential Bayesian inference of the unobserved states via the Kalman Filter. This in turn allows the estimation of the innovations or prediction error which in turn can be used to estimate the likelihood of the model. Once the likelihood has been obtained estimation of model parameters and comparison of models is possible.

Figure 5

Connectivity matrix used for large scale simulations. Sample estimated nerve fiber pathways between voxels of the left and right occipital poles (left panel) of cortical tessellation shown in Figure. From these pathways anatomical connectivity values between these voxels were estimated. For use in the large scale neural mass simulations shown in Figure 2a, in reality all the voxels on both cortical and thalamic surfaces were used to calculate the overall connectivity matrix. This connectivity matrix (dimensions 16138 × 16138) is summarized for purposes of illustration here to a region to regional representation corresponding to 90 well‐known anatomical areas (right panel). The same technique also estimates the length of the fiber connections that is used to infer conduction delays between voxels. Note the extreme sparseness of connections.