Dynamic Granger causality based on Kalman filter for evaluation of functional network connectivity in fMRI data

Abstract

Increasing interest in understanding dynamic interactions of brain neural networks leads to formulation of sophisticated connectivity analysis methods. Recent studies have applied Granger causality based on standard multivariate autoregressive (MAR) modeling to assess the brain connectivity. Nevertheless, one important flaw of this commonly proposed method is that it requires the analyzed time series to be stationary, whereas such assumption is mostly violated due to the weakly nonstationary nature of functional magnetic resonance imaging (fMRI) time series. Therefore, we propose an approach to dynamic Granger causality in the frequency domain for evaluating functional network connectivity in fMRI data. The effectiveness and robustness of the dynamic approach was significantly improved by combining a forward and backward Kalman filter that improved estimates compared to the standard time-invariant MAR modeling. In our method, the functional networks were first detected by independent component analysis (ICA), a computational method for separating a multivariate signal into maximally independent components. Then the measure of Granger causality was evaluated using generalized partial directed coherence that is suitable for bivariate as well as multivariate data. Moreover, this metric provides identification of causal relation in frequency domain, which allows one to distinguish the frequency components related to the experimental paradigm. The procedure of evaluating Granger causality via dynamic MAR was demonstrated on simulated time series as well as on two sets of group fMRI data collected during an auditory sensorimotor (SM) or auditory oddball discrimination (AOD) tasks. Finally, a comparison with the results obtained from a standard time-invariant MAR model was provided.

Figure 1. Simulated data and comparison of dynamic and standard MAR approach

Illustration of non-stationary fMRI time series simulation process, where time series A causes time series B (a); spectrograms of GPDC (both directions) obtained by performing dynamic MAR model based on Kalman filtering before and after smoothing procedure (b); median frequency profiles of GPDC for dynamic MAR and time-invariant MAR estimates (c). The dashed line shows upper bound of 95 % confidence interval under the hypothesis of no connectivity between the components.

Figure 2. Experimental paradigms

Experimental paradigms for SM task (left) and AOD task (right).

Figure 3. ICA spatial maps and associated time courses

Axial slices of selected group spatial maps of independent components, including brain activation of interest, and their associated averaged time courses (over all subjects) for SM (left) and AOD task (right).

Figure 4. Group results of SM task

Set of group results of the dynamic GPDC estimates for SM task. The spectrograms of GPDC, associated with particular causal direction between components, represent the median over all subjects and the single frequency profile is the median over time axis of the associate spectrogram. The dashed line shows upper bound of 95 % confidence interval under the hypothesis of no connectivity between the components. The thin line represents GPDC spectra estimates based on standard MAR model.

Figure 5. Group results of AOD task

Set of group results of the dynamic GPDC estimates for AOD task. The spectrograms of GPDC, associated with particular causal direction between components, represent the median over all subjects and the single frequency profile is the median over time axis of the associate spectrogram. The dashed line shows upper bound of 95 % confidence interval under the hypothesis of no connectivity between the components. The thin line represents GPDC spectra estimates based on standard MAR model.

Figure 6. Connectivity results of functional networks integrations

Directional integrations of functional networks for SM (left) and AOD (right) task, which were revealed as a significant after dynamic MAR estimation, at task-related frequencies 0.031 Hz and 0.033 Hz, respectively. The values above arrows represent magnitudes of GPDC measure at particular driving frequencies. The thin dashed line represents connectivity structure identified by the standard MAR model (without pertinence constrain to single frequency).

Figure 7. Comparison of GPDC based on dynamic and standard MAR as an effect of different model order

The example of GPDC spectra estimates for different model orders (3; 5; 7; 8; 10; 12; 14; 16) obtained by using dynamic MAR (top) and standard time-invariant MAR (bottom), applied to ICA time course of A1 and IFG activations from SM data. Displayed spectra were calculated as a median over all subjects and runs. Black triangle marks correctly estimated spectral peaks that are related to task frequency of fMRI data. White triangle marks wrong identification of causal direction at particular frequency.

Figure 8. Model order selection based on AIC and BIC, and the effect of using ICA time courses compared to original fMRI time courses on estimated GPDC spectra

Estimated BIC and AIC criteria functions of model order for bivariate (2 time courses) and multivariate (4 time courses) standard MAR model applied to SM and AOD data (a). Displayed criteria functions represent the median over all subjects and runs, and their minima, which were used as a suggested model order, are marked. The comparison of resulting GPDC spectra, when the ICA time courses and the time courses from original data (the mean and the 1^st principle component) were used (b).