HHS Public Access Author manuscript Author Manuscript
Adv Appl Stat. Author manuscript; available in PMC 2015 December 07. Published in final edited form as: Adv Appl Stat. 2014 ; 40(1): 61–74.
PARTICLE FILTERING WITH SEQUENTIAL PARAMETER LEARNING FOR NONLINEAR BOLD fMRI SIGNALS Jing Xia1 and Michelle Yongmei Wang2,* 1Public
Health Sciences Division, Fred Hutchinson Cancer Research Center, Seattle, WA, U. S.
A.
Author Manuscript
2Departments
of Statistics, Psychology and Bioengineering, University of Illinois at UrbanaChampaign Champaign, IL, U. S. A.
Abstract
Author Manuscript
Analyzing the blood oxygenation level dependent (BOLD) effect in the functional magnetic resonance imaging (fMRI) is typically based on recent ground-breaking time series analysis techniques. This work represents a significant improvement over existing approaches to system identification using nonlinear hemodynamic models. It is important for three reasons. First, instead of using linearized approximations of the dynamics, we present a nonlinear filtering based on the sequential Monte Carlo method to capture the inherent nonlinearities in the physiological system. Second, we simultaneously estimate the hidden physiological states and the system parameters through particle filtering with sequential parameter learning to fully take advantage of the dynamic information of the BOLD signals. Third, during the unknown static parameter learning, we employ the low-dimensional sufficient statistics for efficiency and avoiding potential degeneration of the parameters. The performance of the proposed method is validated using both the simulated data and real BOLD fMRI data.
Keywords particle filtering; BOLD fMRI; nonlinear dynamics; estimation
1. Introduction Author Manuscript
Functional magnetic resonance imaging (fMRI) has been a standard tool for visualizing regional brain activations during cognitive tasks. The widely used fMRI is based on the blood oxygenation level dependent (BOLD) signal, which can reflect the sluggish, nonlinear and non-local hemodynamic response due to the involved neuronal activity. A standard way of analyzing the BOLD fMRI is the general linear model (GLM) method [1]. In GLM, the brain responses to the stimuli are assumed linear and the whole brain is considered as a black box characterized by a canonical hemodynamic response function (HRF). The GLM approach may be suboptimal for fMRI data analysis for the reason that the brain responses are nonlinear to the functional stimuli. To consider the comprehensive
*
Corresponding author: [email protected].
Xia and Wang
Page 2
Author Manuscript
properties of BOLD fMRI, a dynamic model including all underlying physiological variables and hidden physiological states is highly preferable. The balloon model for the BOLD fMRI proposed by Buxton [2, 3] can describe the nonlinear neural activities in the brain. It is formulated as a mixed continuous discrete time system of nonlinear stochastic differential equations. Riera et al. [4] incorporated previous model and the hemodynamic model which links the stimulus sequence and flow-inducing signal to make a more comprehensive nonlinear model for BOLD signals. These models have a high impact on fMRI analysis; but how to fit and estimate the underlying state-space variables and parameters remains an active area of research.
Author Manuscript Author Manuscript
Non-linear filtering theory has many applications, especially for estimation of continuous stochastic status from non-additive noise. Kalman filtering and extended Kalman filtering have been widely applied in the estimation of stochastic differential equations [5], but they only deal with Gaussian noise and have limitations in approximating the inherent nonlinearities of the dynamics. The linearization approaches [4] also can be applied for simplifying complex nonlinear differential equations, though they cannot provide a sufficient description of the nonlinear nature of a dynamic system, either, as has been demonstrated in [6]. Recently, Johnston et al. applied a particle filtering method for simultaneous estimation of the hidden physiological states and the system parameters based on an iterative coordinate descent approach [6, 7]. However, for simplicity, the involved system parameters are treated as fixed values instead of random variables during estimation, which may not fully take advantage of the complete sequential information as time increases. In addition, different optimization methods and criteria are used for the estimation of the involved parameters, which are also not optimized simultaneously. The maximum likelihood based on which the parameters are updated thus has difficulty in reaching the global maximum. Hettiarachchi et al. used an auxiliary particle filter and smoothing method to jointly estimate BOLD signal hidden physiological states and system parameters [8]. However, the parameter estimation might not be converged because of poor identifiability of the model. The methods used in [9] have similar limitations.
Author Manuscript
In this paper, a more comprehensive particle filtering method with sequential parameter learning is proposed and applied to the BOLD fMRI signals. Particle filtering is a sequential Monte Carlo method based on point mass representations of probability densities. Underlying system parameters are modeled as stationary and unknown, though how to handle the presence of unknown static parameters is still challenging [10, 11] and has yet been applied to the BOLD fMRI. The non-dynamic of the parameters tends to make the parameter samples degenerate into one or a few different values as time increases. Here, we apply a low-dimensional sufficient statistics instead of high-dimensional hidden states for sequential parameter learning, in which a linear computational cost is achieved for efficiency and the degenerate problem is avoided. Please note that in this work, we share the same motivation as Johnston et al.; however, the methods are different. We use Bayesian framework, with all static parameters following some specific distributions and the posterior distributions updated based on measurements. We draw samples that provide the parameter distributions; the sample means are then used to approximate the expected values for the static parameters.
Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 3
Author Manuscript
2. BOLD Signal Modeling The BOLD signal model describes the changes in cerebral oxygen extraction fraction and cerebral blood volume during brain activation. The hemodynamic balloon model [3] links the relationship among the BOLD signal, yt, and the intrinsic state variables including the flow-inducing signal, st, the cerebral blood flow (CBF), ft, the cerebral blood volume (CBV), vt, and the total deoxyhemoglobin (Hb), qt. In this model, the venous compartment is considered to be a balloon-like structure, inflated by increased blood flow resulting from an increase in local neuronal activity. The blood flow increases more than the cerebral metabolic rate of oxygen, thus locally reducing the deoxyhemoglobin content which in turn causes a slight increase in the local MR signal, creating the BOLD effect. The measured BOLD signal yt is modeled as a nonlinear function of ‘snapshots’ of the continuously evolving states with additive white Gaussian noise et:
Author Manuscript
(1)
where V0 is resting blood volume fraction, and {k1 = 7E0, k2 = 2, k3 = 2E0 − 0.3} are generally accepted coefficients, with E0 the resting oxygen extraction fraction. The dynamical model consists of four non-linear differential equations: (2)
(3)
Author Manuscript
(4)
(5)
where ut is a task stimulus; ε is the neuronal efficacy; τs is the signal decay time constant; τf is the autoregulatory time constant; τ0 is the mean transit time of venous compartment at rest; and α is the stiffness component.
Author Manuscript
Equations (1) to (5) describe the relation among input ut, output measurement yt and a set of hidden state variables xt = [st, ft, vt, qt]T. The above model was extended in [4] to include the noise in the dynamic system: (6)
where wt is a Wiener process. F(·) encompasses all the right-hand sides of the 4 differential equations. Note that in this work, we also model the variance for each of the hidden state
Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 4
Author Manuscript
variables and dwt is then a vector . This is not achieved in previous work such as [4, 6, 7]. Thus, the complete set of parameters for the BOLD signal to be estimated .
is:
There is no analytic solution for the system of nonlinear ordinary differential equations, by applying the Euler-Maruyama method [12] with sampling at integer multiples of Δt and based on equation (6), the discrete-time dynamics can be written as: (7)
3. Particle Filtering with Sequential Parameter Learning Author Manuscript
Consider a general state-space model with the measurement yt, the unobserved hidden state variable xt, and a parameter vector θ. Let Δt = 1, the model is specified in terms of the probability densities:
Author Manuscript
The initial state distribution p(x0|θ) is assumed to be known. In the BOLD model, the joint filtering and parameter distribution is unavailable in closed form and shall be approximated using Monte Carlo methods. The required posterior density function can be represented by a set of random samples with associated weights, and all estimators are then computed by these samples. When the parameters θ are assumed to be known, we have the pure filtering problem. In this case, the goal is to generate samples {xt} from the filtering density p(xt | y1:t) for each time t based on importance sampling and resampling [5, 7, 11]. The posterior filtered density p(xt| y1:t) can be approximated as: (8)
It can be shown that as sample size N → ∞, the approximation approaches the true posterior density p(xt|y1:t). If the prior density p(xt|xt−1) is chosen as the importance, the weight is
Author Manuscript
given by
.
We propose to use the following particle filtering strategy to deal with the presence of unknown parameters. When a new observation yt arrives, the particles are updated in order to represent the new posterior p(x1:t|y1:t). In many cases, only p(xt | y1:t) is of interest. In order to handle the presence of unknown static parameters, the usual trick is to include the parameters as part of the state vector {xt, θ}. Because of the non-dynamic feature of the parameters, samples of θ at time t can only take the values given at time t − 1. However, Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 5
Author Manuscript
most values become very unlikely when new observations are available; this will result in an impoverishment of the set of distinct θ values. The distribution of θ depends on x1:t, y1:t, whose dimension gets higher as t increases. Thus we use some sufficient statistics Tt = Tt(x1:t, y1:t), a low dimension function, to update the posterior distribution of θ. The use of sufficient statistics can make the filter more efficient as time evolves (the details are provided in the next two paragraphs). The particle filtering is also used in [6, 7] for BOLD fMRI; however, the parameter updating is based on the maximum likelihood estimation (MLE) using knowledge of the updated state sequence estimates, which is a simplified approach and may not be optimal for dynamic systems. We employ the sequential parameter learning algorithm to better utilize the dynamic system information for improved robustness, which is summarized in Table 1.
Author Manuscript
Assume the priors for all variances as: (σf)2 ~ IG(υ, φf1), (σs)2 ~ IG(υ, φs1), (σv)2 ~ IG(υ, φv1), (σq)2 ~ IG(υ, φq1), and (σe)2 ~ IG(υ, φe1), where IG is the inverse Gamma distribution. The posterior distributions of the variance are thus:
The sufficient statistics for υt, φft, and so on are updated according to:
Author Manuscript
Regarding the static parameters, ε; τs, τf, τ0, E0, V0, we set them as generalized beta distributions; the posterior distributions for these parameters are still generalized beta. For example, τt ~ Beta(0, 5, αt, βt), E0 = Beta(0.2, 0.8, αt, βEt) and V0 = Beta(0.01, 0.04, αt, βVt). The updated sufficient statistics are αt = α0 + t − 1, βt = β0 + 2υv1 + t, βEt = βE0 + 2υq1 + t, βVt = βV0 + 2υe1 + t · υv1, where υq1 and υe1 are hyper parameters of the priors for the variance of vt, qt and et.
4. Experiments Author Manuscript
4.1. Simulated data In this section, some examples of the dynamic model will be examined in order to evaluate the performance of the presented particle filtering when static parameters are known or unknown. Ground truth against which to test the state estimation algorithm is not attainable through real fMRI data, for which reason we simulate the BOLD signal model and generate the block design data as in [6]. The stimulus sequence is used to drive the physiological variables’ responses. Simulation of hidden status requires continuous time state dynamics;
Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 6
Author Manuscript
we apply the Euler-Maruyama method. When Δt → 0, the discrete time dynamics converge in mean square sense to the continuous time dynamics. The discrete time dynamics is: xt+Δt = xt + F(xt, ut, θ) · Δt + wt+Δt. The output BOLD signal could be generated by using equation (1) and
, which accounts for measurement and instrumentation noise.
4.1.1. Comparison of extended Kalman filter (EKF) and particle filtering— Extended Kalman filter (EKF) is a local linearization of the dynamic equations but the model parameters need to be known when using the EKF. The hemodynamic model parameters in the simulated BOLD data were chosen as ε = 1, τs = 1.25, τf = 2.5, E0 = 0.5, V0 = 0.025 and the sampling interval Δt = 1. The variances for all the state functions are,
Author Manuscript
, , , , and τ0 = 2.5. Figure respectively: 1 demonstrates that the particle filter maintains accurate state estimates while the EKF fails to track the state sequences precisely. In order to quantify the performance, the standard mean squared error (MSE) is used as the evaluation criterion, as in [4]. In the block design, the MSEs are, respectively, 3.5546 (for EKF) and 0.9210 (for particle filter). The eventrelated design data we tested have the slimier performance as the block design. 4.1.2. Comparison of non-filtered fMRI and particle filtered signals—In this experiment, we would like to test the robustness and sensitivity of the proposed particle filtering method on the simulated block design data with different signal-to-noise ratios (SNRs). Since the BOLD signal and the CBF are simultaneously measurable in physiological experiments [2], here we compare the GLM-based t-map results from filtered CBF (ft) and from non-filtered original fMRI data (yt). The estimation for system parameters is also provided. The empirical variance of the noise
and that of the signal
Author Manuscript
that various SNRs (SNR = −2, −1, 0; experiments.
are set so
) are obtained for
Figure 2 shows the estimation results of the hidden states and the parameters. The activation maps in Figure 2(a) demonstrate that our particle filtering method is not sensitive to the noise, leading to consistent activation maps despite the variations in SNR. However, the GLM results directly based on the original fMRI measurement get worse as SNR decreases. The estimated patterns for the hidden states are significantly different between the activation and non-activation regions (Figure 2(c) shows all the estimated states for the activation voxels but only two most important states (CBF ft and CBV vt) for the non-activation voxels due to the space limit). Also, the estimation of the system parameters converges within reasonable iterations (Figure 2(b)).
Author Manuscript
4.2. Real fMRI data The real fMRI data (single subject, size: 53*63*46*360) is obtained from the SPM data site (http://www.fil.ion.ucl.ac.uk/~wpenny/datasets/attention.html) with a visual motion task. The subject was scanned during four runs with 90 image volumes in each run. Four conditions - ‘fixation’, ‘attention’, ‘no attention’ and ‘stationary’- are used. Figure 3 shows particle filtering results of the hidden states and some estimated system parameters, for the corresponding voxels in V1 Right, V5 Right, PP Right (right posterior parietal cortex) and Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 7
Author Manuscript
PFC Right (right dorsolateral prefrontal cortex) regions. As can be seen, the filtered hidden states for the four regions share similar activation patterns. Also, the baseline of CBF (ft) in PFC Right is lower than that in V1 Right, implying the activation is stronger in V1 Right [3]. These results are consistent with the interpretations in [14]. The resting parameters, t0, E0 and V0 converge well and share similar values over different regions.
5. Conclusions
Author Manuscript
In this paper, we presented a comprehensive particle filtering method with sequential parameter learning and applied the proposed approach to the nonlinear BOLD fMRI signal analysis. Superior performance of the method over the extended Kalman filter is demonstrated; the successful joint state and parameter estimation is also illustrated in the experiments. Our proposed particle filtering provides an alternative way of analyzing the nonlinear and stochastic brain hemodynamic system, and in fact is ideally suited for handling the nonlinearity inherent in the BOLD signals and the propagating uncertainty associated with the sequential Monte Carlo method. Future directions of our work include: incorporating the spatial information into the dynamic system analysis, comparing our method with the one in [6] based on the iterative coordinate descent technique, and further validation of the proposed methods.
Acknowledgement This work was supported in part by NIH grant K25AG033725.
References Author Manuscript Author Manuscript
[1]. Friston KJ, Holmes AP, Worsley KJ, Poline JP, Frith CD, Frackowiak RSJ. Statistical parametric maps in functional imaging: a general linear approach. Human Brain Mapping. 1995; 2:189–210. [2]. Buxton RB, Wong EC, Frank LR. Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magnetic Resonance in Medicine. 1998; 39(6):855–864. [PubMed: 9621908] [3]. Buxton RB, Uludag K, Dubowitz DJ, Liu TT. Modeling the hemodynamic response to brain activation. NeuroImage. 2004; 23:S220–S233. [PubMed: 15501093] [4]. Riera JJ, Watanabe J, Kazuki I, Naoki M, Aubert E, Ozaki T, Kawashima R. A state-space model of the hemodynamic approach: nonlinear filtering of BOLD signals. NeuroImage. 2004; 21:547– 567. [PubMed: 14980557] [5]. Arulampalam MS, Maskell S, Gordon N, Clapp T. A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Trans. Signal Processing. 2002; 50:174–188. [6]. Johnston LA, Duff E, Mareels I, Egan GF. Nonlinear estimation of the BOLD signal. NeuroImage. 2008; 40:504–514. [PubMed: 18203623] [7]. Johnston, LA.; Duff, E.; Egan, GF. Particle filtering for nonlinear BOLD signal analysis; Int. Conf. Medical Image Computing and Computer Assisted Intervention; 2006; p. 292-299. [8]. Hettiarachchi, IT.; Mohamed, S.; Nahavandi, S. Identification of nonlinear fMRI models using auxiliary particle filter and kernel smoothing method; Annual International Conference of the IEEE Engineering in Medicine and Biology Society; 2012; p. 4212-4216. [9]. Murry, L.; Storkey, A. Advances in Neural Information Processing Systems. 2007. Continuous time particle filtering for fMRI. [10]. Polson NG, Stroud JR, Muller P. Practical filtering with sequential parameter learning. J. Royal Statistical Society Series B. 2008; 70:413–428.
Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 8
Author Manuscript
[11]. Storvik G. Particle filters in state space models with the presence of unknown static parameters. IEEE Trans. Signal Processing. 2002; 50:281–289. [12]. Kloeden, PE.; Platen, E. Numerical Solution of Stochastic Differential Equations. Springer; 1999. [13]. Friston KJ, Mechelli A, Turner R, Price CJ. Nonlinear responses in fMRI: the balloon model, Volterra kernels, and other hemodynamics. NeuroImage. 2000; 12(4):466–477. [PubMed: 10988040] [14]. Buchel C, Friston KJ. Modulation of connectivity in visual pathways by attention: cortical inferences evaluated with structural equation modeling and fMRI. Cerebral Cortex. 1997; 7:768– 778. [PubMed: 9408041]
Author Manuscript Author Manuscript Author Manuscript Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 9
Author Manuscript Figure 1.
Author Manuscript
Performance comparison of EKF and our particle filter for the simulated block design data. (Solid line: ground truth; thick dashed line: results from the particle filter; dotted line: results from EKF. When not visible, the particle filtering coincides with the ground truth.)
Author Manuscript Author Manuscript Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 10
Author Manuscript Author Manuscript
Figure 2.
Performance of the present particle filtering for the simulated block design data: (a) comparison of activation maps for various SNRs (top: from fMRI measurement; bottom: from the filtered CBF), (b) estimation of the system parameters for the activation voxel when SNR = −1 (horizontal line: ground truth), (c) estimation of the BOLD signal and the hidden states for activation (with all estimated states shown) and non-activation (with only CBF and CBV shown here) voxels when SNR = −1 (solid line: ground truth; thick dashed line: estimated values).
Author Manuscript Author Manuscript Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 11
Author Manuscript Figure 3.
Author Manuscript
Particle filtering results for the corresponding voxels in V1 Right and PFC Right regions (two of the four runs are shown), (a) positions of the respective regions, (b) the stimulus signal (top), the measurement of the BOLD signal (bottom, dotted line), and the estimated BOLD signal (bottom, solid line), (c) estimation of the hidden states, (d) estimation of some parameters relating to the resting state.
Author Manuscript Author Manuscript Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 12
Table 1
Author Manuscript
The proposed particle filter algorithm with sequential parameter learning at each time t For the samples i = 1, …, N ----Sample �~f
(� ∣ Tit−1, yt ) f (� ∣ Tit−1, yt ) is proposal distributions for
θ, which can be updated with time.
i
----Sample xt ~p
(xt ∣ xit−1). i
----Calculate the importance weights: wt
(
i
= p yt ∣ xt
)
End For
N
Calculate total weight: total
= SUM
{wti }i=1s
For i = 1, …, Ns
Author Manuscript
i
----Normalize: wt
= total −1wti
End For Resample algorithm
N {xit , wti }t=1
Author Manuscript Author Manuscript Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Adv Appl Stat. Author manuscript; available in PMC 2015 December 07. Published in final edited form as: Adv Appl Stat. 2014 ; 40(1): 61–74.
PARTICLE FILTERING WITH SEQUENTIAL PARAMETER LEARNING FOR NONLINEAR BOLD fMRI SIGNALS Jing Xia1 and Michelle Yongmei Wang2,* 1Public
Health Sciences Division, Fred Hutchinson Cancer Research Center, Seattle, WA, U. S.
A.
Author Manuscript
2Departments
of Statistics, Psychology and Bioengineering, University of Illinois at UrbanaChampaign Champaign, IL, U. S. A.
Abstract
Author Manuscript
Analyzing the blood oxygenation level dependent (BOLD) effect in the functional magnetic resonance imaging (fMRI) is typically based on recent ground-breaking time series analysis techniques. This work represents a significant improvement over existing approaches to system identification using nonlinear hemodynamic models. It is important for three reasons. First, instead of using linearized approximations of the dynamics, we present a nonlinear filtering based on the sequential Monte Carlo method to capture the inherent nonlinearities in the physiological system. Second, we simultaneously estimate the hidden physiological states and the system parameters through particle filtering with sequential parameter learning to fully take advantage of the dynamic information of the BOLD signals. Third, during the unknown static parameter learning, we employ the low-dimensional sufficient statistics for efficiency and avoiding potential degeneration of the parameters. The performance of the proposed method is validated using both the simulated data and real BOLD fMRI data.
Keywords particle filtering; BOLD fMRI; nonlinear dynamics; estimation
1. Introduction Author Manuscript
Functional magnetic resonance imaging (fMRI) has been a standard tool for visualizing regional brain activations during cognitive tasks. The widely used fMRI is based on the blood oxygenation level dependent (BOLD) signal, which can reflect the sluggish, nonlinear and non-local hemodynamic response due to the involved neuronal activity. A standard way of analyzing the BOLD fMRI is the general linear model (GLM) method [1]. In GLM, the brain responses to the stimuli are assumed linear and the whole brain is considered as a black box characterized by a canonical hemodynamic response function (HRF). The GLM approach may be suboptimal for fMRI data analysis for the reason that the brain responses are nonlinear to the functional stimuli. To consider the comprehensive
*
Corresponding author: [email protected].
Xia and Wang
Page 2
Author Manuscript
properties of BOLD fMRI, a dynamic model including all underlying physiological variables and hidden physiological states is highly preferable. The balloon model for the BOLD fMRI proposed by Buxton [2, 3] can describe the nonlinear neural activities in the brain. It is formulated as a mixed continuous discrete time system of nonlinear stochastic differential equations. Riera et al. [4] incorporated previous model and the hemodynamic model which links the stimulus sequence and flow-inducing signal to make a more comprehensive nonlinear model for BOLD signals. These models have a high impact on fMRI analysis; but how to fit and estimate the underlying state-space variables and parameters remains an active area of research.
Author Manuscript Author Manuscript
Non-linear filtering theory has many applications, especially for estimation of continuous stochastic status from non-additive noise. Kalman filtering and extended Kalman filtering have been widely applied in the estimation of stochastic differential equations [5], but they only deal with Gaussian noise and have limitations in approximating the inherent nonlinearities of the dynamics. The linearization approaches [4] also can be applied for simplifying complex nonlinear differential equations, though they cannot provide a sufficient description of the nonlinear nature of a dynamic system, either, as has been demonstrated in [6]. Recently, Johnston et al. applied a particle filtering method for simultaneous estimation of the hidden physiological states and the system parameters based on an iterative coordinate descent approach [6, 7]. However, for simplicity, the involved system parameters are treated as fixed values instead of random variables during estimation, which may not fully take advantage of the complete sequential information as time increases. In addition, different optimization methods and criteria are used for the estimation of the involved parameters, which are also not optimized simultaneously. The maximum likelihood based on which the parameters are updated thus has difficulty in reaching the global maximum. Hettiarachchi et al. used an auxiliary particle filter and smoothing method to jointly estimate BOLD signal hidden physiological states and system parameters [8]. However, the parameter estimation might not be converged because of poor identifiability of the model. The methods used in [9] have similar limitations.
Author Manuscript
In this paper, a more comprehensive particle filtering method with sequential parameter learning is proposed and applied to the BOLD fMRI signals. Particle filtering is a sequential Monte Carlo method based on point mass representations of probability densities. Underlying system parameters are modeled as stationary and unknown, though how to handle the presence of unknown static parameters is still challenging [10, 11] and has yet been applied to the BOLD fMRI. The non-dynamic of the parameters tends to make the parameter samples degenerate into one or a few different values as time increases. Here, we apply a low-dimensional sufficient statistics instead of high-dimensional hidden states for sequential parameter learning, in which a linear computational cost is achieved for efficiency and the degenerate problem is avoided. Please note that in this work, we share the same motivation as Johnston et al.; however, the methods are different. We use Bayesian framework, with all static parameters following some specific distributions and the posterior distributions updated based on measurements. We draw samples that provide the parameter distributions; the sample means are then used to approximate the expected values for the static parameters.
Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 3
Author Manuscript
2. BOLD Signal Modeling The BOLD signal model describes the changes in cerebral oxygen extraction fraction and cerebral blood volume during brain activation. The hemodynamic balloon model [3] links the relationship among the BOLD signal, yt, and the intrinsic state variables including the flow-inducing signal, st, the cerebral blood flow (CBF), ft, the cerebral blood volume (CBV), vt, and the total deoxyhemoglobin (Hb), qt. In this model, the venous compartment is considered to be a balloon-like structure, inflated by increased blood flow resulting from an increase in local neuronal activity. The blood flow increases more than the cerebral metabolic rate of oxygen, thus locally reducing the deoxyhemoglobin content which in turn causes a slight increase in the local MR signal, creating the BOLD effect. The measured BOLD signal yt is modeled as a nonlinear function of ‘snapshots’ of the continuously evolving states with additive white Gaussian noise et:
Author Manuscript
(1)
where V0 is resting blood volume fraction, and {k1 = 7E0, k2 = 2, k3 = 2E0 − 0.3} are generally accepted coefficients, with E0 the resting oxygen extraction fraction. The dynamical model consists of four non-linear differential equations: (2)
(3)
Author Manuscript
(4)
(5)
where ut is a task stimulus; ε is the neuronal efficacy; τs is the signal decay time constant; τf is the autoregulatory time constant; τ0 is the mean transit time of venous compartment at rest; and α is the stiffness component.
Author Manuscript
Equations (1) to (5) describe the relation among input ut, output measurement yt and a set of hidden state variables xt = [st, ft, vt, qt]T. The above model was extended in [4] to include the noise in the dynamic system: (6)
where wt is a Wiener process. F(·) encompasses all the right-hand sides of the 4 differential equations. Note that in this work, we also model the variance for each of the hidden state
Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 4
Author Manuscript
variables and dwt is then a vector . This is not achieved in previous work such as [4, 6, 7]. Thus, the complete set of parameters for the BOLD signal to be estimated .
is:
There is no analytic solution for the system of nonlinear ordinary differential equations, by applying the Euler-Maruyama method [12] with sampling at integer multiples of Δt and based on equation (6), the discrete-time dynamics can be written as: (7)
3. Particle Filtering with Sequential Parameter Learning Author Manuscript
Consider a general state-space model with the measurement yt, the unobserved hidden state variable xt, and a parameter vector θ. Let Δt = 1, the model is specified in terms of the probability densities:
Author Manuscript
The initial state distribution p(x0|θ) is assumed to be known. In the BOLD model, the joint filtering and parameter distribution is unavailable in closed form and shall be approximated using Monte Carlo methods. The required posterior density function can be represented by a set of random samples with associated weights, and all estimators are then computed by these samples. When the parameters θ are assumed to be known, we have the pure filtering problem. In this case, the goal is to generate samples {xt} from the filtering density p(xt | y1:t) for each time t based on importance sampling and resampling [5, 7, 11]. The posterior filtered density p(xt| y1:t) can be approximated as: (8)
It can be shown that as sample size N → ∞, the approximation approaches the true posterior density p(xt|y1:t). If the prior density p(xt|xt−1) is chosen as the importance, the weight is
Author Manuscript
given by
.
We propose to use the following particle filtering strategy to deal with the presence of unknown parameters. When a new observation yt arrives, the particles are updated in order to represent the new posterior p(x1:t|y1:t). In many cases, only p(xt | y1:t) is of interest. In order to handle the presence of unknown static parameters, the usual trick is to include the parameters as part of the state vector {xt, θ}. Because of the non-dynamic feature of the parameters, samples of θ at time t can only take the values given at time t − 1. However, Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 5
Author Manuscript
most values become very unlikely when new observations are available; this will result in an impoverishment of the set of distinct θ values. The distribution of θ depends on x1:t, y1:t, whose dimension gets higher as t increases. Thus we use some sufficient statistics Tt = Tt(x1:t, y1:t), a low dimension function, to update the posterior distribution of θ. The use of sufficient statistics can make the filter more efficient as time evolves (the details are provided in the next two paragraphs). The particle filtering is also used in [6, 7] for BOLD fMRI; however, the parameter updating is based on the maximum likelihood estimation (MLE) using knowledge of the updated state sequence estimates, which is a simplified approach and may not be optimal for dynamic systems. We employ the sequential parameter learning algorithm to better utilize the dynamic system information for improved robustness, which is summarized in Table 1.
Author Manuscript
Assume the priors for all variances as: (σf)2 ~ IG(υ, φf1), (σs)2 ~ IG(υ, φs1), (σv)2 ~ IG(υ, φv1), (σq)2 ~ IG(υ, φq1), and (σe)2 ~ IG(υ, φe1), where IG is the inverse Gamma distribution. The posterior distributions of the variance are thus:
The sufficient statistics for υt, φft, and so on are updated according to:
Author Manuscript
Regarding the static parameters, ε; τs, τf, τ0, E0, V0, we set them as generalized beta distributions; the posterior distributions for these parameters are still generalized beta. For example, τt ~ Beta(0, 5, αt, βt), E0 = Beta(0.2, 0.8, αt, βEt) and V0 = Beta(0.01, 0.04, αt, βVt). The updated sufficient statistics are αt = α0 + t − 1, βt = β0 + 2υv1 + t, βEt = βE0 + 2υq1 + t, βVt = βV0 + 2υe1 + t · υv1, where υq1 and υe1 are hyper parameters of the priors for the variance of vt, qt and et.
4. Experiments Author Manuscript
4.1. Simulated data In this section, some examples of the dynamic model will be examined in order to evaluate the performance of the presented particle filtering when static parameters are known or unknown. Ground truth against which to test the state estimation algorithm is not attainable through real fMRI data, for which reason we simulate the BOLD signal model and generate the block design data as in [6]. The stimulus sequence is used to drive the physiological variables’ responses. Simulation of hidden status requires continuous time state dynamics;
Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 6
Author Manuscript
we apply the Euler-Maruyama method. When Δt → 0, the discrete time dynamics converge in mean square sense to the continuous time dynamics. The discrete time dynamics is: xt+Δt = xt + F(xt, ut, θ) · Δt + wt+Δt. The output BOLD signal could be generated by using equation (1) and
, which accounts for measurement and instrumentation noise.
4.1.1. Comparison of extended Kalman filter (EKF) and particle filtering— Extended Kalman filter (EKF) is a local linearization of the dynamic equations but the model parameters need to be known when using the EKF. The hemodynamic model parameters in the simulated BOLD data were chosen as ε = 1, τs = 1.25, τf = 2.5, E0 = 0.5, V0 = 0.025 and the sampling interval Δt = 1. The variances for all the state functions are,
Author Manuscript
, , , , and τ0 = 2.5. Figure respectively: 1 demonstrates that the particle filter maintains accurate state estimates while the EKF fails to track the state sequences precisely. In order to quantify the performance, the standard mean squared error (MSE) is used as the evaluation criterion, as in [4]. In the block design, the MSEs are, respectively, 3.5546 (for EKF) and 0.9210 (for particle filter). The eventrelated design data we tested have the slimier performance as the block design. 4.1.2. Comparison of non-filtered fMRI and particle filtered signals—In this experiment, we would like to test the robustness and sensitivity of the proposed particle filtering method on the simulated block design data with different signal-to-noise ratios (SNRs). Since the BOLD signal and the CBF are simultaneously measurable in physiological experiments [2], here we compare the GLM-based t-map results from filtered CBF (ft) and from non-filtered original fMRI data (yt). The estimation for system parameters is also provided. The empirical variance of the noise
and that of the signal
Author Manuscript
that various SNRs (SNR = −2, −1, 0; experiments.
are set so
) are obtained for
Figure 2 shows the estimation results of the hidden states and the parameters. The activation maps in Figure 2(a) demonstrate that our particle filtering method is not sensitive to the noise, leading to consistent activation maps despite the variations in SNR. However, the GLM results directly based on the original fMRI measurement get worse as SNR decreases. The estimated patterns for the hidden states are significantly different between the activation and non-activation regions (Figure 2(c) shows all the estimated states for the activation voxels but only two most important states (CBF ft and CBV vt) for the non-activation voxels due to the space limit). Also, the estimation of the system parameters converges within reasonable iterations (Figure 2(b)).
Author Manuscript
4.2. Real fMRI data The real fMRI data (single subject, size: 53*63*46*360) is obtained from the SPM data site (http://www.fil.ion.ucl.ac.uk/~wpenny/datasets/attention.html) with a visual motion task. The subject was scanned during four runs with 90 image volumes in each run. Four conditions - ‘fixation’, ‘attention’, ‘no attention’ and ‘stationary’- are used. Figure 3 shows particle filtering results of the hidden states and some estimated system parameters, for the corresponding voxels in V1 Right, V5 Right, PP Right (right posterior parietal cortex) and Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 7
Author Manuscript
PFC Right (right dorsolateral prefrontal cortex) regions. As can be seen, the filtered hidden states for the four regions share similar activation patterns. Also, the baseline of CBF (ft) in PFC Right is lower than that in V1 Right, implying the activation is stronger in V1 Right [3]. These results are consistent with the interpretations in [14]. The resting parameters, t0, E0 and V0 converge well and share similar values over different regions.
5. Conclusions
Author Manuscript
In this paper, we presented a comprehensive particle filtering method with sequential parameter learning and applied the proposed approach to the nonlinear BOLD fMRI signal analysis. Superior performance of the method over the extended Kalman filter is demonstrated; the successful joint state and parameter estimation is also illustrated in the experiments. Our proposed particle filtering provides an alternative way of analyzing the nonlinear and stochastic brain hemodynamic system, and in fact is ideally suited for handling the nonlinearity inherent in the BOLD signals and the propagating uncertainty associated with the sequential Monte Carlo method. Future directions of our work include: incorporating the spatial information into the dynamic system analysis, comparing our method with the one in [6] based on the iterative coordinate descent technique, and further validation of the proposed methods.
Acknowledgement This work was supported in part by NIH grant K25AG033725.
References Author Manuscript Author Manuscript
[1]. Friston KJ, Holmes AP, Worsley KJ, Poline JP, Frith CD, Frackowiak RSJ. Statistical parametric maps in functional imaging: a general linear approach. Human Brain Mapping. 1995; 2:189–210. [2]. Buxton RB, Wong EC, Frank LR. Dynamics of blood flow and oxygenation changes during brain activation: the balloon model. Magnetic Resonance in Medicine. 1998; 39(6):855–864. [PubMed: 9621908] [3]. Buxton RB, Uludag K, Dubowitz DJ, Liu TT. Modeling the hemodynamic response to brain activation. NeuroImage. 2004; 23:S220–S233. [PubMed: 15501093] [4]. Riera JJ, Watanabe J, Kazuki I, Naoki M, Aubert E, Ozaki T, Kawashima R. A state-space model of the hemodynamic approach: nonlinear filtering of BOLD signals. NeuroImage. 2004; 21:547– 567. [PubMed: 14980557] [5]. Arulampalam MS, Maskell S, Gordon N, Clapp T. A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Trans. Signal Processing. 2002; 50:174–188. [6]. Johnston LA, Duff E, Mareels I, Egan GF. Nonlinear estimation of the BOLD signal. NeuroImage. 2008; 40:504–514. [PubMed: 18203623] [7]. Johnston, LA.; Duff, E.; Egan, GF. Particle filtering for nonlinear BOLD signal analysis; Int. Conf. Medical Image Computing and Computer Assisted Intervention; 2006; p. 292-299. [8]. Hettiarachchi, IT.; Mohamed, S.; Nahavandi, S. Identification of nonlinear fMRI models using auxiliary particle filter and kernel smoothing method; Annual International Conference of the IEEE Engineering in Medicine and Biology Society; 2012; p. 4212-4216. [9]. Murry, L.; Storkey, A. Advances in Neural Information Processing Systems. 2007. Continuous time particle filtering for fMRI. [10]. Polson NG, Stroud JR, Muller P. Practical filtering with sequential parameter learning. J. Royal Statistical Society Series B. 2008; 70:413–428.
Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 8
Author Manuscript
[11]. Storvik G. Particle filters in state space models with the presence of unknown static parameters. IEEE Trans. Signal Processing. 2002; 50:281–289. [12]. Kloeden, PE.; Platen, E. Numerical Solution of Stochastic Differential Equations. Springer; 1999. [13]. Friston KJ, Mechelli A, Turner R, Price CJ. Nonlinear responses in fMRI: the balloon model, Volterra kernels, and other hemodynamics. NeuroImage. 2000; 12(4):466–477. [PubMed: 10988040] [14]. Buchel C, Friston KJ. Modulation of connectivity in visual pathways by attention: cortical inferences evaluated with structural equation modeling and fMRI. Cerebral Cortex. 1997; 7:768– 778. [PubMed: 9408041]
Author Manuscript Author Manuscript Author Manuscript Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 9
Author Manuscript Figure 1.
Author Manuscript
Performance comparison of EKF and our particle filter for the simulated block design data. (Solid line: ground truth; thick dashed line: results from the particle filter; dotted line: results from EKF. When not visible, the particle filtering coincides with the ground truth.)
Author Manuscript Author Manuscript Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 10
Author Manuscript Author Manuscript
Figure 2.
Performance of the present particle filtering for the simulated block design data: (a) comparison of activation maps for various SNRs (top: from fMRI measurement; bottom: from the filtered CBF), (b) estimation of the system parameters for the activation voxel when SNR = −1 (horizontal line: ground truth), (c) estimation of the BOLD signal and the hidden states for activation (with all estimated states shown) and non-activation (with only CBF and CBV shown here) voxels when SNR = −1 (solid line: ground truth; thick dashed line: estimated values).
Author Manuscript Author Manuscript Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 11
Author Manuscript Figure 3.
Author Manuscript
Particle filtering results for the corresponding voxels in V1 Right and PFC Right regions (two of the four runs are shown), (a) positions of the respective regions, (b) the stimulus signal (top), the measurement of the BOLD signal (bottom, dotted line), and the estimated BOLD signal (bottom, solid line), (c) estimation of the hidden states, (d) estimation of some parameters relating to the resting state.
Author Manuscript Author Manuscript Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
Xia and Wang
Page 12
Table 1
Author Manuscript
The proposed particle filter algorithm with sequential parameter learning at each time t For the samples i = 1, …, N ----Sample �~f
(� ∣ Tit−1, yt ) f (� ∣ Tit−1, yt ) is proposal distributions for
θ, which can be updated with time.
i
----Sample xt ~p
(xt ∣ xit−1). i
----Calculate the importance weights: wt
(
i
= p yt ∣ xt
)
End For
N
Calculate total weight: total
= SUM
{wti }i=1s
For i = 1, …, Ns
Author Manuscript
i
----Normalize: wt
= total −1wti
End For Resample algorithm
N {xit , wti }t=1
Author Manuscript Author Manuscript Adv Appl Stat. Author manuscript; available in PMC 2015 December 07.
