Synopsis

The signal fluctuations in Functional Magnetic Resonance Imaging (fMRI) have been proved suitable for investigating brain connectivity. Sequential Monte Carlo methods (“particle filters”) aim at estimating internal states in dynamic systems when only partial and noisy observations are available. Differently from the majority of techniques commonly used for investigating brain connectivity, which assume stationarity, particle filters are designed for stochastic time-varying systems. We present a particle filtering algorithm tailored to fMRI data, whose purpose is to help assessing the presence of causal influences that certain brain areas may exert over others. The algorithm has been validated on a simulated network and applied to real fMRI data.

Introduction

Theory, Materials and Methods

A time-varying first-order auto-regressive model was adopted to represent the relationship across fMRI timecourses (“observation equation”, Fig.1/Eq.1).

The time-varying nature of the model is accounted for by the a_ri,t coefficients, that represent the causal influence exerted by area i over area r. The a_ri,t coefficients are assumed to evolve in time according to a first-order Markov-chain model (“state equation”, Fig.1/Eq.2).

The output of the algorithm are the coefficients a_ri,t, which are predicted via the Bayes rule when observations become available (Fig.1/Eq.3). The posterior distribution is approximated with a finite set of N weighted samples (“particles”) by using a recently published particle filter¹, adapted for fMRI timecourses. In this study, N was set to 2000.

To validate the algorithm, synthetic data of length=100 and six nodes following the model in Fig.1/Eq.1&2 were generated with MATLAB, both without additive Gaussian noise and with SNR=6dB. The time-averaged a_ri coefficients estimated by particle filtering were compared to the actual inputs of the simulated data, and their relationship was assessed by Pearson’s correlation.

Two actual fMRI datasets, previously acquired on two healthy volunteers with a GE MR950 7T scanner, were used to test the proposed method. Data were acquired with TR=2s and consisted of 240 volumes with isotropic voxels of size (1.5mm)³. During acquisition, the subjects’ thumb- and index-fingertips were stimulated in a block-design fashion (20s stimulation ON, 20s stimulation OFF) by a tactile stimulation device (Linari Engineering). During the stimulation period, subjects were asked to move the stimulated fingers. Networks of four nodes were studied. Nodes consisted of the voxels in Regions-of-Interest (ROIs) covering primary somatosensory (S1), primary motor (M1), supplementary motor and parietal cortex. ROIs were manually drawn on each subject on one slice only, to avoid potential slice timing confounds. The optimal order of the autoregressive model describing the timeseries was 1, as estimated by the Schwartz criterion². In these fMRI datasets, the a_ri coefficients estimated by particle filtering were compared to the delayed correlation c_ri between signals x_r,t and x_i,t-1, which reflect the time-invariant causal influence exerted by node i over node r. a_ri coefficients were considered significantly different from zero (p<0.05) when their modulus was higher than a threshold value=0.1 determined on the basis of a null distribution of coefficients.

Results and Discussion

Fig.2 demonstrates a significant correlation (both p<0.001) between the a_ri coefficients estimated by particle filtering and the actual inputs of the simulated noiseless and noisy data (Pearson’s correlation coefficient=0.96 and 0.59, respectively).

Fig.3 shows that the non-zero a_ri coefficients estimated by particle filtering in actual fMRI data correlated with the time-delayed correlation coefficients c_ri. Pearson’s correlation coefficient was 0.41.

Across time, only a few a_ri,t (including those representing the causal influences of M1 over itself and supplementary motor area over parietal cortex) exhibited, in both subjects, statistically significant difference between frames with stimulus/task and those during rest (t-test p<0.05). The temporal oscillations of the coefficient representing the causal influence of S1 over M1 passed the statistical threshold only in one subject (p<0.05) and not in the other (p=0.080). Factors that might have prevented the detection of further stimulus-locked oscillations (if they were actually present) include time resolution and/or the length of the task blocks.

A similar approach has previously been applied to resting state fMRI³, however the matrix of a_ri coefficients was assumed symmetric (a_ri=a_ir). Further validation will assess in which scenarios the proposed approach might be used to capture transient changes in effective connectivity, both in resting state and in stimulus/task-driven fMRI.

Acknowledgements

References

1. Ancherbak, S., Kuruoglu, E. E. & Vingron, M. Time-Dependent Gene Network Modelling by Sequential Monte Carlo. IEEE/ACM Trans Comput Biol Bioinform 13, 1183–1193 (2016).

2. Roebroeck, A., Formisano, E. & Goebel, R. Mapping directed influence over the brain using Granger causality and fMRI. NeuroImage 25, 230–242 (2005).

3. Ahmad, M. F., Murphy, J., Vatansever, D., Stamatakis, E. A. & Godsill, S. Tracking changes in functional connectivity of brain networks from resting-state fMRI using particle filters. in 798–802 (IEEE, 2015). doi:10.1109/ICASSP.2015.7178079