Real-Time fMRI with Turbo-BrainVoyager 2.0
I feel very grateful when looking back to the events of this year’s summer: It started with the OHBM conference in Florence and the BrainVoyager 10th anniversary party where the winners of the first BrainVoyager Design Awards were announced. Then there was the exciting Football World Cup in Germany, which made especially my Italian friends very happy :-). This was followed by a wonderful vacation in Sri Lanka and Maldives and a trip to the Formula 1 Grand Prix at Monza. And there were many nice short trips, e.g. to Oxford (ASSC10), Warsaw (ESMRMB 2006) and Chieti (ISBET 2006), Despite so many exciting events, there was still some time for writing and programming. Here I report about my work on Turbo-BrainVoyager.
Turbo-BrainVoyager, or “TBV” for short, allows to analyze fMRI data in real-time, i.e. during an ongoing experiment. This movie provides an impression about how TBV is working. While analyzing the data, TBV remains highly interactive allowing the user to explore interesting functional brain regions as they emerge during real-time analysis. Turbo-BrainVoyager 2.0 adds several new features since the last release (v1.8) making it suitable for advanced applications such as neurofeedback. One of the advanced new features is the possibility to visualize dynamic statistical maps on 3D anatomical data sets as well as on cortex meshes. These representations may even reside in AC-PC or Talairach space. In BrainVoyager QX, the whole 4D time course data (FMR-STC) is transformed to AC-PC or Talairach space (VMR-VTC data) before calculating and visualizing contrast maps or ROIs. This normalization step is especially useful for efficient integration of data across subjects. Since TBV analyzes single subject data only, AC-PC or Talairach normalization is performed in a different way. Statistical analyses are always performed in the original slice-based (FMR-STC) space. To render statistical data on the right positions in AC-PC or Talairach space, each voxel or vertex retrieves the data from the slice-based space using a single spatial transformation step describing how a point in slice space maps to a point in AC-PC or Talairach space and vice versa. When the user, for example, clicks on a cortex mesh in Talairach space, the coordinates of the hit vertex are retrieved in Talairach 3D space and then sent “backward” through the spatial transformation matrix indexing the corresponding voxel in slice space. A little problem is given by the fact that the indexed slice coordinates will not be of integral value but will index coordinate values falling in-between the original voxel grid. This problem is typically solved by interpolation using a weighted sum of neighboring voxels to estimate the value at the indexed position. In TBV, trilinear interpolation is used since it provides good results and because it can be implemented efficiently.
Another area of important improvements concerns statistical analysis tools including dynamic display of estimated beta values for all active ROIs and the possibility to turn on/off contrast maps during real-time processing. It is now also possible to inspect and modify ROIs at the voxel level using a zoomed slice visualization, which is very useful for neurofeedback applications. Besides 3D motion correction, it is now also possible to apply a Gaussian spatial smoothing step to the raw data. While these preprocessing steps take extra time, the amount of time is constant for each newly arriving volume allowing to guarantee real-time processing for arbitrarily long runs. This is, unfortunately, not the case for standard General Linear Model analyses, which must be performed for each voxel using the whole time course data to estimate optimal beta values. A GLM can be written as:
where y refers to the time course of a voxel, X refers to the design matrix and e refers to the errors (residuals). The unknown vector b is determined using the standard least squares approach minimizing the sum of squared error values. The respective beta values can be calculated as follows:
This equation requires matrix multiplication of the design matrix with its transpose followed by matrix inversion and the result is multiplied with the result of the matrix-vector product of the design matrix and the data. With longer and longer time courses in a real-time situation, the time to calculate the betas will steadily increase over time when using this equation because the design matrix and data vector increases steadily. One trick to speed up computation is to calculate the first term (matrix inversion of X’X) only once since this part does not change from voxel to voxel. While this approach, used in BrainVoyager QX, reduces overall computation time substantially, the time required to perform statistics will still increase somewhat with each new data point. Even with fast computer hardware, there will be a time point reached when the calculations per volume will require more time than the time available, which is the time of one TR (typically 1 - 3 seconds) minus the time for other calculations. At this point real-time performance will be lost and the analysis will lack behind the data. Fortunately this problem can be avoided by using a version of the GLM, which allows calculating beta values incrementally. Instead of using a standard least squares procedure to find the right beta values, TBV uses a recursive least squares (RLS) approach to update the betas from time point to time point using only the newly arriving data. The equation to do that looks at first more complex then the one of the standard GLM:
In this recursive equation, the new beta values (bt+1) are calculated by starting with the last beta values adding a term, which can be calculated by using only the new values of the design matrix at that time point (xt+1). There is, however, still the (X’X)-1 term. Fortunately this term can be also updated recursively using only the information of the new values (single “row” xt+1) of the design matrix:
Like in the case of calculating the new betas, the calculation of the new inverted matrix uses the last matrix to calculate the new one using only the new data of the design matrix. Since these recursive equations use a fixed amount of calculation at each time point, real-time performance is guaranteed for arbitrary long time series. For completeness, here is the equation for the f term, which is calculated first since it is used in both update equations:
Note that the update equation for the (X’X)-1 term does again not depend on the data and can thus be calculated once at the beginning of each time point. Thus only the beta update equation has to be calculated for each voxel. This strategy is used in TBV to make GLM analyses very fast and to guarantee a constant amount of computation time from time point to time point.