Incremental activation detection for real-time fMRI series using robust kalman filter.
Functional magnetic resonance imaging (fMRI) offers a noninvasive method in studying human brain functions by recording blood-oxygen-level-dependent (BOLD) signal changes related to neuronal activity across the brain with high spatial resolution . Real-time fMRI (rt-fMRI) is a method to assess the acquired data for evidence of an experimentally induced effect at every intracerebra voxel individually and simultaneously. In rt-fMRI, data are processed as fast as they are acquired . For real-time fMRI applications, mapping the activations within a repetition time makes it possible to interact with fMRI experiments in a much more efficient way . Online functional mapping enables researchers to monitor data quality, evolve experimental protocols more rapidly, perform interactive experimental paradigms for neurological investigation , achieve neurofeedback by providing feedback ofbrain activation to the subject in real time 5], which may have potential use in clinical applications .
In common fMRI experiments, MRI scanner acquires whole brain data at an interval of 2 seconds, also called repetition time. To meet the real-time requirements, all the processing steps of real-time fMRI need to be completed within a repetition time. Simple real-time fMRI processing steps consist of data reconstruction, spatial realignment (head motion correction), and statistical analysis. Among them, incremental statistical analysis on each voxel of the fMRI dataset will result in huge computational costs. To overcome the computational costs of the statistical analysis, a number of incremental activation detection algorithms have been developed for rt-fMRI applications.
Cox et al.  proposed the first real-time incremental activation detection algorithm and the correlation based activation detection methods are not able to model multiple experimental and confounding effects simultaneously. Gembris et al.  proposed correlation analysis method using sliding window technique. Friston et al.  proposed the general linear model (GLM), which can be used as a unified framework in the analysis of fMRI data and support multiple experimental design, but it cannot be used to rt-fMRI applications, because it needs all of the data to do the statistical analysis. The widely used correlation based techniques are special cases of GLM with a white noise model for the temporal
errors in the signal. Bagarinao et al.  proposed an algorithm using an orthogonalization procedure to estimate the coefficient of general linear models. Roche et al.  proposed an algorithm using extended kalman filter (EKF) method to fit a general linear model on fMRI time courses. This technique adopts the GLM-AR model under the assumption that the fMRI noise is significantly autocorrelated. The extended kalman filter is able to fit incrementally GLM coefficients along with the one-order autoregressive noise model. In this paper, we mainly focus on the EKF method, because it requires low computation costs and memory costs, and the design matrix can be assembled incrementally, which makes it possible for more complex interactive experiment design.
The fMRI time series has low signal to noise ratio . The noise in the fMRI signal is complicated. It includes not only the usual MRI sensor noise but also physiological fluctuations that affect the signal. The fluctuations of MRI scanner or severe head motion may generate sparse noise in the signal. In our model, under the assumption that the sparse noise will not change the noise distribution, a variation is added to the extended kalman filter to handle the additional measurement noise term, that is, sparse; this term can be used to model the sparse measurement outliers.
In recent years, with the developments of convex optimization, Mattingley and Boyd  created a robust kalman filter by replacing the standard measurement update, which can be interpreted as the result of solving a similar convex minimization problem, which includes an [l.sub.1] term to handle the sparse noise. We use the robust kalman filter method to solve the value of the sparse term in our model.
Real-time fMRI signal is three-dimensional volume data at each scan during a repetition time and the intensity of each voxelrepresentsthe blood oxygen level associatedwithneural activities. Data of each incoming fMRI scans is spatially aligned with the first scan of the series. Voxel values at different scan time points are arranged to time sequence, forming the measured time series. Each voxel will form a time series, and the length of time series is growing with the time increasing. The time series of each voxel can be calculated independently, so in the following discussion we only consider the situation of a single voxel time series.
2.1. Extended Kalman Filter Incremental Detection. Roche et al.  proposed an extended kalman filter based algorithm to fit incrementally the general linear model along with a one-order autoregressive noise model.
The general linear model (GLM) explains the measured time series y = [[[y.sub.1], [y.sub.2], ..., [y.sub.n]].sup.T] in terms of a linear combination of the explanatory variables [X.sub.p] = [[x.sub.1p], [x.sub.2p], ..., [x.sub.np]].sup.T] plus an error term. The explanatory variable contains paradigm-related regressors based on the experiment design and signal characteristics and regressors are obtained by convolving the different stimulation onsets with a canonical hemodynamic response function. Model the low-frequency drift, hence enabling us to "detrend" the signal (we use polynomials up to order three). Then combine the explanatory variables into a design matrix X = [[X.sub.1], [X.sub.2], ..., [X.sub.p]]. The GLM can be expressed as
y = X[beta] + [epsilon]. (1)
The design matrix is X = [[X.sub.1], [X.sub.2], ..., [X.sub.p]], where [X.sub.p] is the explanatory variable, [beta] = [[[b.sub.1], [b.sub.2], ..., [b.sub.p]].sup.T] contains (unknown) parameters which represent the coefficients of the explanatory variables. The design matrix contains paradigm related regressors based on the experiment design. The regressors can be obtained by convolving the different stimulation onsets with a canonical hemodynamic response function, or modeling the low-frequency drift to detrend the signal.
In the GLM-AR model, it is assumed that e is a stationary Gaussian zero mean AR(1) random process and the relationship between [[epsilon].sub.i] and [[epsilon].sub.i-1] can be expressed as.
[[epsilon].sub.i] = a[[epsilon].sub.i-1] + [n.sub.i]. (2)
a is an (unknown) autocorrelation parameter, and [n.sub.i] is a white noise with instantaneous Gaussian distribution N(0, [[sigma].sup.2]).
Alexis Roche proved that the maximum likelihood estimate of ([beta], a) can be computed independently from a and they found that one can reach the maximum of the likelihood function by finding the minimum of [[rho].sub.i]([beta], a):
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3)
where [x.sub.i] = [[[x.sub.i1], [x.sub.i2], ..., [x.sub.ip]].sup.T] denote the values of explanatory variables at time i.
First, they combine the parameters to be estimated [beta] and a into a (p + 1) x 1 state vector b = [[[beta]; a].sup.T] and assume that, at time i, the current estimate of b is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].
Secondly, linearize the error function [[rho].sub.i]([beta], a) to [[rho].sub.i](b) around the current estimate using a first-order Taylor expansion:
min [[rho].sub.i] (b) = [q.sub.i] - [u.sup.T.sub.i]b, (4)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].
Finally, they solved the (nonlinear) least-squares regression problem by means of an EKF. The EKF updates the parameters using the following recursion:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (5)
where [k.sub.i] is the kalman gain in ith step and [[SIGMA].sub.i] denote the normalized posterior covariance matrix of b given the information available at time i.
At time i, the posterior covariance matrix of the state vector is approximated by Cov([b.sub.i]) [approximately equal to] [[??].sup.2.sub.i][[SIGMA].sub.i] and the incremental update rule for a is as follows:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (6)
The number of explanatory variables is p, and [beta] is the p x 1 vector which contains the coefficients related to the explanatory variables. [[SIGMA].sup.pxP.sub.i] is obtained by extracting the left superior p x p block from the matrix [[??].sup.2.sub.i][[SIGMA].sub.i]. For a given contrast vector c, we can identify the voxels that show a contrasted effect:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (7)
where [t.sub.i] = [([c.sup.T][[SIGMA].sup.pxP.sub.i]c).sup.-1/2][c.sup.T] p and erf (*) is the error function. Testing for positive activations can be achieved at any time i by thresholding the image of [t.sub.i]-statistics.
2.2. Outlier Detection Method. The sensor failures or measurement outliers will cause the sparse measurement noise and they may cause rapid degrade on the detection performance. We derive a new algorithm to detect the outliers in order to eliminate the effect on the kalman filter algorithm.
We suppose that there is a sparse term [z.sub.i] which is always zero, and it has no effect on the distribution of the total noise; then the general linear model can be modified as
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (8)
Then linearize the [[rho].sub.i]([beta], a, [z.sub.i]) to [[rho].sub.i](b, [z.sub.i]): [rho](b, [z.sub.i]) [approximately equal to] [q.sub.i] - [u.sup.T.sub.i] b - [z.sub.i], (9)
where [q.sub.i] = [y.sub.i] - [a.sub.i-1] [x.sup.T.sub.i-1] [[beta].sub.i-1] and [u.sub.i] = [[[x.sub.i] - [a.sub.i-1] [x.sup.T.sub.i-1], [y.sub.i-1] - [x.sup.T.sub.i-1][[beta].sub.i-1]-[z.sub.i-1]].sup.T], [z.sub.i] is an unknown variable at time i, while [z.sub.i-1] is an known variable at time i.
The linearized constraint equation is approximate as follows:
[q.sub.i] = [u.sup.T.sub.i][b.sub.i] + [z.sub.i] + [v.sub.i]. (10)
The measurement update step of standard kalman filter algorithm is essentially an optimization problem, and the linearized parameter optimization problem can be described as follow:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)
where b and [v.sub.i] are the unknown parameters to be estimated, Z denotes the steady-state error covariance associated with predicting the next state b, and the measurement is [q.sub.i] and the measurement noise term [v.sub.i] is a white noise with instantaneous Gaussian distribution N(0, V).
We use the robust kalman filter method to detect the outlier hidden in the measurement by replacing the measurement update with the solution of a similar convex minimization problem, which includes an [l.sub.1] term to handle the sparse noise.
To (approximately) handle the additional sparse noise term [z.sub.i], we modify the kalman filter measurement update step. The modified optimization problem is as follows:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (12)
In the optimization problem (12), b, [z.sub.i], and [v.sub.i] are the variables to be estimated.
To solve this convex optimization problem, we adopt a fast transform method proposed by Mattingely and Boyd 13], which makes problem solving become more effective:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (13)
After the transformation, the original problem was transformed into an equivalent convex quadratic program problem:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (14)
With variable [z.sub.i], solving the convex quadratic program, we can achieve the sparse noise [z.sub.i] at each time. Fortunately, for one voxel at time i, the size of [e.sub.i], [z.sub.i] and [Q.sub.i] is 1 x 1. Hence this optimization problem is equivalent to solve the solution of a piecewise-quadratic function. This problem has an analytical solution, means we can use analytical expression instead of the searching optimization loop.
2.3. Robust Extended Kalman Filter. The standard kalman filter consists of alternating time and measurement updates. Since [beta] is slowly varying parameter, so it remains unchanged in the time update step. Both of the time and measurement updates are derived by the minimum mean-square error estimates of b. As is shown in (5), the measurement update equation is as follows:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (15)
After solving the problem, we achieve an estimate outlier value [[??].sub.i] and then combine it into the update step:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (16)
Furthermore, we combine the [??] into the incremental update rule for [sigma]:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (17)
For the real-time i-test, the method is exactly the same as Alexis Roche's method.
Finally, the algorithm recursion is summarized as follows:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (18)
The algorithm is tested on a single run from an fMRI experiment involving a visual and auditory task. The protocol is block design, and the run consisted of 10 blocks, each block including one activation epoch (20 s) and one control epoch (10 s). We aim at finding the voxels associated with the visual and auditory function; each activation epoch, the visual stimulus, and auditory stimuli are present, but the intensity is different. The data has an abrupt head motion during the scan in the 45th repetition time (TR) and the head motion caused severe motion artifacts. The repetition time was 2 seconds for a total of 152 scans. Functional images have 80 * 80 * 33 voxels.
As is shown in Figure 1, the red curve is the reference vector obtained by convolving the stimulation onset with a canonical hemodynamic response function, which is assumed the time series should be; the blue curve is the time series data of an active voxel; obviously there is an outlier at the 45th TR.
In Figure 2, all of the three curves tend to be stabilized and the correlation coefficient [beta] of activate voxel is all converge to around 26. The green curve presents Cox method, and it has a great drop at the outlier scan, and the extended kalman filter method has an abnormal rise at the 45th TR, while our method has no great fluctuation at the outlier scan. This shows its robustness to the sparse noise.
Cox method derived a threshold for active detection, so there is no T value in their method. Figure 3 shows the comparison of the i-test value between the Robust extended kalman filter and the extended kalman filter. We can see that both of the two algorithms have significant decrease on the outlier. Obviously, our method is more stable at the outlier scan and has a higher T value after the outlier scan.
We also tested the algorithm on an inactivate voxel. In Figure 4, the blue curve is the time series of an inactive voxel and at the 45th TR there is also an outlier.
The t-test value of the inactive voxel should be low, as is shown in Figure 5; the correlation curves of the three methods tend to be stable below zero. There is a big fluctuation at the outlier scan according to Cox method and our method finally obtained a stable value lower than the extended kalman method. Before the outlier scan, the correlations of the two methods show little difference. After the outlier scan, our method becomes lower than the extended kalman method.
In Figure 6, a comparison of the T value between our method and the extended kalman filter method is shown. We can see that T value of both the two algorithms has fluctuations at the outlier, and the fluctuations of our method are smaller than those of extended kalman filter method. Our method tends to be stable at a lower level. This again shows our new method's robustness to the sparse noise in the inactive voxel.
Figure 7 is the result of outlier detection, as is shown in Figures 1 and 4; the original signal contains a big outlier at the 45th TR and the margin of the two outlier is different. Our outlier detection algorithm detect two outliers from both the active and inactive voxel at 45th TR, it is shown to be robust for both the active and inactive voxels.
We have applied our detection algorithm to the whole brain data, in the case where threshold value the same, the final activation maps are shown in Figure 8; left half is the activation result achieved by the extended kalman algorithm and right half is the result of the algorithm described in this paper. As expected for a visual auditory task, activations are found in the visual cortex and auditory cortex on both sides. This means that, under the same threshold of T value, we can get more voxels associated with the task.
We have a C++ version of the algorithm, and we test it on a 16-core processor workstation, in which computing hole brain (80 *80* 32) voxels costs 0.3 s to 0.4 s; there is plenty of time to do some other processes. Our algorithm can be used on the real-time fMRI application.
In our algorithm, parameter [lambda] needs to be defined before the experiment. In the optimization object function, the parameter [lambda] can be regarded as the weight of the sparse term. With a large [lambda], the outlier detection algorithm will find nothing and result zero in [z.sub.i], the measurement update will be the same as the extended kalman filter method. So if no outlier is detected, the result of the algorithm will be the same as that of in the extended kalman filter. Small [lambda] will make the algorithm modify the measurement update frequently. Under this situation, the mistaken outlier will have an effect on the noise distribution and cause the unstability of the kalman filter. Here we can regard [lambda] as the detection threshold and by holding a large threshold we can detect the obvious outlier, which may be caused by the machine failure. The value of [lambda] in our algorithm is in a wide range. We suggest taking a relatively large [lambda] to test if the result is stable. In this experiment, we test [lambda] from 50 to 300 and the algorithm can detect the outlier and the result is robust to the outlier.
The fMRI signal has a low signal to noise ratio, so we cannot give an exact threshold value and this algorithm provides a method to threshold the signal, and we can even achieve the value of the outlier and use it to modify the detection algorithm. fMRI experiment is complex, which needs cooperation with subject and operator; any problem of them may cause the experiment to fail. Our algorithm can detect this kind of failure and it can detect the outliers, which may have a great effect on follow processes, where we can stop the experiment to check the problem or just mark it and eliminate the outlier data after the experiment.
In future works, we will explore a proper method to determine X, just enough to detect the outlier but not cause the unstability of the kalman filter algorithm.
The robust kalman filter is introduced in the activation detection of fMRI experiment. Convex optimization method is used to modify the extended kalman filter by introducing a sparse noise term. The robustness of our method to the sparse noise is improved. Moreover, the performance of the proposed method does not degrade rapidly when disturbances are involved. When applied to the time series voxels, our method can obtain more stable i-test value in both activate and inactivate voxels. The proposed algorithm can also detect outliers, which may have a great effect on following processes. The detected outlier information can be used to reject or modify the bad data, which may have potential benefit for real-time applications that require higher data quality.
This work is supported by the National High Technology Research and Development Program of China (no. 2012AA011603).
 N. K. Logothetis, "MR imaging in the non-human primate: studies of function and of dynamic connectivity," Current Opinion in Neurobiology, vol. 13, no. 5, pp. 630-642, 2003.
 N. Weiskopf, R. Sitaram, O. Josephs et al., "Real-time functional magnetic resonance imaging: methods and applications," Magnetic Resonance Imaging, vol. 25, no. 6, pp. 989-1003, 2007
 R. Christopher deCharms, "Applications of real-time fMRI," Nature Reviews Neuroscience, vol. 9, no. 9, pp. 720-729, 2008.
 S. M. LaConte, "Decoding fMRI brain states in real-time," NeuroImage, vol. 56, no. 2, pp. 440-454, 2011.
 N. Weiskopf, "Real-time fMRI and its application to neurofeedback," NeuroImage, vol. 62, no. 2, pp. 682-692, 2011.
 N. Birbaumer, A. R. Murguialday, C. Weber, and P. Montoya, "Chapter 8 neurofeedback and brain-computer interface: clinical applications," International Review of Neurobiology, vol. 86, pp. 107-117, 2009.
 R. W. Cox, A. Jesmanowicz, and J. S. Hyde, "Real-time functional magnetic resonance imaging," Magnetic Resonance in Medicine, vol. 33, no. 2, pp. 230-236, 1995.
 D. Gembris, J. G. Taylor, S. Schor, W. Frings, D. Suter, and S. Posse, "Functional magnetic resonance imaging in real time (FIRE): sliding-window correlation analysis and referencevector optimization," Magnetic Resonance in Medicine, vol. 43, no. 2, pp. 259-268, 2000.
 K. J. Friston, A. P. Holmes, K. J. Worsley, J.-P Poline, C. D. Frith, and R. S. J. Frackowiak, "Statistical parametric maps in functional imaging: a general linear approach," Human Brain Mapping, vol. 2, no. 4, pp. 189-210, 1994.
 E. Bagarinao, K. Matsuo, T Nakai, and S. Sato, "Estimation of general linear model coefficients for real-time application," NeuroImage, vol. 19, no. 2, pp. 422-429, 2003.
 A. Roche, P J. Lahaye, and J. B. Poline, "Incremental activation detection in fMRI series using Kalman filtering," in Proceedings of the IEEE International Symposium on Biomedical Imaging: Nano to Macro, vol. 1, pp. 376-379, April 2004.
 B. B. Averbeck and D. Lee, "Effects of noise correlations on information encoding and decoding," Journal of Neurophysiology, vol. 95, no. 6, pp. 3633-3644, 2006.
 J. Mattingely and S. Boyd, "Real-time convex optimization in signal processing," IEEE Signal Processing Magazine, vol. 27, no. 3, pp. 50-61, 2010.
Liang Li, Bin Yan, Li Tong, Linyuan Wang, and Jianxin Li
China National Digital Switching System Engineering and Technological Research Center, Zheng Zhou 450002, China
Correspondence should be addressed to Li Tong; firstname.lastname@example.org
Received 27 September 2013; Accepted 6 December 2013; Published 6 January 2014
Academic Editor: Sabri Arik
|Printer friendly Cite/link Email Feedback|
|Title Annotation:||Research Article|
|Author:||Li, Liang; Yan, Bin; Tong, Li; Wang, Linyuan; Li, Jianxin|
|Publication:||Computational and Mathematical Methods in Medicine|
|Date:||Jan 1, 2014|
|Previous Article:||A mixture modeling framework for differential analysis of high-throughput data.|
|Next Article:||A harmonic linear dynamical system for prominent ECG feature extraction.|