Computing moment-to-moment BOLD activation for real-time neurofeedback

Abstract

Estimating moment-to-moment changes in blood oxygenation level dependent (BOLD) activation levels from functional magnetic resonance imaging (fMRI) data has applications for learned regulation of regional activation, brain state monitoring, and brain-machine interfaces. In each of these contexts, accurate estimation of the BOLD signal in as little time as possible is desired. This is a challenging problem due to the low signal-to-noise ratio of fMRI data. Previous methods for real-time fMRI analysis have either sacrificed the ability to compute moment-to-moment activation changes by averaging several acquisitions into a single activation estimate or have sacrificed accuracy by failing to account for prominent sources of noise in the fMRI signal. Here we present a new method for computing the amount of activation present in a single fMRI acquisition that separates moment-to-moment changes in the fMRI signal intensity attributable to neural sources from those due to noise, resulting in a feedback signal more reflective of neural activation. This method computes an incremental general linear model fit to the fMRI time series, which is used to calculate the expected signal intensity at each new acquisition. The difference between the measured intensity and the expected intensity is scaled by the variance of the estimator in order to transform this residual difference into a statistic. Both synthetic and real data were used to validate this method and compare it to the only other published real-time fMRI method.

Figure 1

Schematic demonstrating the moment to moment activation estimation method. The activation estimate at time t is the difference between the data at t and the model reconstruction using only nuisance bases However, the full model fit $\widehat{y}$ is still computed 1) so that the residual error in the model fit can be used to estimate the efficiency of the voxel signal and 2) so that the timecourse variance associated with task-related neural signals is not attributed to drift, which would bias the model fit.

Figure 2

Example of synthetic fMRI timecourse generation. The top three panels show individual components of the resulting signal. The top panel shows the contribution to the final signal from hemodynamic response to activation (1% signal change). The second panel shows the contribution from scanner drift or physiological noise, measured from real resting fMRI scans (in this case 0.5% signal change). The third panel shows the Gaussian white noise vector generated for this example (SNR of 2). The large panel shows the resulting fMRI timecourse data (the sum of the first three panels and a constant 500 unit offset), the post-hoc GLM reconstruction of this timecourse, and the incremental GLM reconstruction from a fit performed only on data up to that time. Below this, the model difference panel shows the reconstruction error, which is the squared difference between the post-hoc and incremental reconstructions. The bottom panel shows the evolution of the model parameter estimate of the neural basis, which stabilizes near the correct value of 4 after the first stimulus block.

Figure 3

Effect of manipulating SNR of the synthetic fMRI timeseries on the error in model reconstruction between post-hoc and incremental GLM fits. Error bars indicate standard error over the 1000 synthetic timeseries.

Figure 4

Effect of manipulating the strength of white noise and linear drift in the synthetic fMRI timeseries on the error in model reconstruction between post-hoc and incremental GLM fits. The strength of linear drift in % signal change is indicated by bar color, as shown in the legend. Error bars indicate standard error over the 1000 synthetic timeseries.

Figure 5

Effect of manipulating the strength of non-linear drift in the synthetic fMRI timeseries. Bar color indicates % signal change of the drift. Error bars indicate standard error over the mean error from each of the 6 different drift signals.

Figure 6

Effect on GLM reconstruction error of adding temporal autocorrelation to the synthetic fMRI data. The strength of non-linear drift is indicated by bar color using the same color and error bar scheme as in Figure 4 and Figure 5.

Figure 7

Comparison of ground truth neurofeedback with feedback computed using our method with either mean, median, or weighted average ROI combination as well as neurofeedback computed using APSC. The mean t-statistic shown here is the mean regression statistic over all subjects and functional runs of the self-regulation experiment. Methods were compared using a paired t-test, and all differences are significant to p < 0.0001. Error bars represent standard error from the mean.