Hansen, Lars Kai

Jacobsen, Daniel Jakup

Madsen, Kristoffer Hougaard

Crossref

Identification of non-linear models of neural activity in BOLD fMRI

Non-linear hemodynamic models express the BOLD signal as a nonlinear, parametric functional of the temporal sequence of local neural activity. Several models have been proposed for this neural activity. We identify one such parametric model by estimating the distribution of its parameters. These distributions are themselves stochastic, therefore we estimate their variance by epoch based leave-one-out cross validation, using a Metropolis-Hastings algorithm for sampling of the posterior parameter distribution

Online Research Database In Technology

                       General rights Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.  • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain • You may freely distribute the URL identifying the publication in the public portal   If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.     Downloaded from orbit.dtu.dk on: Dec 17, 2017Identification of non-linear models of neural activity in bold fmriJacobsen, Danjal Jakup; Madsen, Kristoffer Hougaard; Hansen, Lars KaiPublished in:3rd IEEE International Symposium on Biomedical Imaging: Macro to NanoLink to article, DOI:10.1109/ISBI.2006.1625077Publication date:2006Document VersionPublisher's PDF, also known as Version of recordLink back to DTU OrbitCitation (APA):Jacobsen, D. J., Madsen, K. H., & Hansen, L. K. (2006). Identification of non-linear models of neural activity inbold fmri. In 3rd IEEE International Symposium on Biomedical Imaging: Macro to Nano (pp. 952-955). IEEE.DOI: 10.1109/ISBI.2006.1625077IDENTIFICATION OF NON-LINEAR MODELS OF NEURAL ACTIVITY IN BOLD FMRIDaniel J. Jacobsen, Kristoffer Hougaard Madsen, Lars Kai HansenIntelligent Signal ProcessingTechnical University of DenmarkABSTRACTNon-linear hemodynamic models express the BOLD signal asa nonlinear, parametric functional of the temporal sequence oflocal neural activity. Several models have been proposed forthis neural activity. We identify one such parametric modelby estimating the distribution of its parameters. These distri-butions are themselves stochastic, therefore we estimate theirvariance by epoch based leave-one-out cross validation, usinga Metropolis-Hastings algorithm for sampling of the posteriorparameter distribution.1. INTRODUCTIONNeuroimaging has made major contributions to our under-standing of the relation between behavior and distribution ofbrain information processing. The richest neuroimaging modal-ity is fMRI based on the so-called BOLD effect, involving thehemodynamic response which is rather sluggish, non-linear,and non-local. With the long term goal of increasing thespatio-temporal resolution of BOLD fMRI, we are interestedin models linking subject behavior, neural activity, the so-called hemodynamic response, and fMRI BOLD observations,see e.g., [1, 2, 3, 4].We will examine the model proposed by Friston et al.[2] which consists of a set of ordinary differential equations(ODE’s) that model the evolution in time of four basic phys-iological state variables: The blood volume v(t), blood in-flow f(t), amount of de-oxyhemoglobine q(t) and a so-called’flow inducing signal’ s(t), collected in the state vector, x(t) =[v(t) q(t) f(t) s(t)]T . The flow inducing signal is driven byan underlying neural activation function ν(t) - a time functiondescribing the local neural activity.The measured BOLD signal yn is then modeled as a non-linear function of ‘snapshots’ of the continuously evolvingstates, with additive white Gaussian noise wn; subscript in-dices are used for these variables to emphasize the discrete‘sampled’ nature.∂x∂t= f(x(t), ν(t))yn = g(x(tn)) + wn(1)The BOLD signal is measured with a sampling interval de-noted TR. The model has seven parameters: σ2w, the vari-ance of wn, and six physiological parameters1, combined inθ =[α & τ0 τs τf E0 σ2w]T . E0 is the so-called ‘resting netoxygen extraction fraction in the capillary bed’.In addition, we assume that the states evolve from an ini-tial known resting state x0 = [0 1 1 1]T . The latter assump-tion is reasonable if a suitably long resting period precedesstimulation sequences. The dynamical model thus consists ofthe set of non-linear differential equations∂v(t)∂t=1τ0(f(t)− v(t)1/α)∂q(t)∂t=1τ0[f(t)1− (1− E0)1/f(t)E0− v(t)(1−α)/αqt]∂s(t)∂t= &ν(t)− s(t)/τs − (f(t)− 1)/τf∂f(t)∂t= s(t)Finally, the BOLD observation model involves the non-linearity,g(x(t)) = V0[(k1 + k2)(1− q(t))− (k2 + k3(1− v(t)))] (2)with a set of empirical constants taking values V0 = 0.02,v0 = 40.3, TE = 0.03, r0 = 25, & = 1.43,k1 = 4.3E0v0TE, k2 = E0&r0TE, k3 = &− 1for BOLD imaging at 1.5T [5]. The BOLD fMRI measure-ments are spatially sampled in volume elements (voxels).Experiments are typically ivided temporally in quasi-independent baseline-activation stimulus ‘epochs’.2. STATISTICAL MODELINGFor given hemodynamic parameters and neural activity, thelikelihood of an epoch is straightforward to set up. First,the hidden states will evolve deterministically according to(1) driven by the given neural activity ν(t). We use a vari-able step-size 4th/5th-order embedded Runge-Kutta methodto solve these [6], with the starting condition x(t = 0) =1we assume the resting blood volume fraction, V0, to be constant at 0.02[2]9520­7803­9577­8/06/$20.00 ©2006 IEEE ISBI 2006Authorized licensed use limited to: D anmarks Tekniske Informationscenter. D ownloaded on N ovember 28 , 2009 at 16:10  from IE E E  X plore .  R estrictions apply. x0, the initial (relaxed) state, x0 = [1 1 1 0]T (all valuesrelative to resting state). This gives a sequence of states,x1:N ! {x1,x2, . . . ,xN}, corresponding to the samplingtimes {t1, t1 + TR, . . . , t1 + N · TR}, where t1 is the start-ing time of the epoch. The mean BOLD signal is given byyn = g(xn; θ)with the observed output modelled asy∗n = yn + wn(θ) (3)As the residuals are assumed normal i.i.d., y∗n ∼ N (yn,σ2w),the likelihood becomesp(y∗1:N |θ) =N∏n=0p(y∗n|θ) =N∏n=0Ny∗n(yn,σ2w).The neural activity function ν(t) is traditionally assumed tobe a simple square wave signal representing signal ‘on’ duringstimulus and signal ‘off’ during baseline [2]. Buxton et al. [3]have recently introduced an alternative dynamical model rep-resenting neural activity which we will use here. This modelpositsν(t) = a(t)− I(t)dIdt=κν(t)− I(t)τ1,(4)where a(t) is the square wave stimulus reference function andI(t) is an inhibitory feedback signal. The values of the param-eters are unknown a priori, although in [3], ranges are given asτ1 ∈ [1; 3], κ ∈ [0; 3]. We refer to them jointly as φ = [κ, τi].Note that the square wave model obtains as a special caseof this non-linear model as κ/τi → 0.In this report we will apply this model to real fMRI dataand investigate the posterior distribution of κ and τi. By wayof leave-one-out cross-validation, the uncertainty on thesedistributions will further be described.3. GENERALIZATION AND ESTIMATIONPROCEDURESThe target for this investigation is the posterior distribution ofthe parameters of the non-linear neural model, φ:p(φ|D) = p(D|φ)p(φ)p(D)=p(D|φ)p(φ)∫p(D|φ)p(φ)dφi.e. the distribution of the parameters conditioned on the ob-served BOLD data, D. This distribution cannot be obtainedanalytically; instead, we employ a Metropolis-HastingsMarkov-chainMonte Carlo (MCMC)method to generate sam-ples from the posterior.3.1. Markov-chain Monte Carlo samplingWe use a Metropolis-Hastings (MH) algorithm [7] startingat an arbitrary state, φ0, and at each step proposing smallchanges in φ from a proposal distribution, in our case a Gaus-sian centered on the ‘current’ state p(φn+1|φn) = N (φn;Σ).The parameters of the hemodynamic model are sampled si-multaneously, but as they are not the focus of this report, theyare not discussed further and ignored for clarity (also the ap-proximate marginal distribution of κ and τi is just the sam-pled values of these, so no further work is required to obtainthe desired distribution).The MHmethod produces samples from the true posteriordistribution in the limit of large number of samples. We em-ploy a set of heuristics to ensure convergence before averag-ing, e.g., inspection of saturation of the training set likelihoodand stabilization of the actual parameter values.A good proposal distribution is major determinant for suc-ces of the Metropolis-Hastings algorithm. Again we invokeheuristics: Starting with a spherical normal distribution of di-mension dimφ, we perform several (short) scout samplingruns. After each of these, the covariance of the generatedsamples is used to adapt the covariance of the proposal, Σ,scaled to give an acceptance rate of around 0.3. This proce-dure greatly speeds up the final sampling run (the samples ofthe initial runs are not used in averaging).For most of the parameters, we use simple uniform pri-ors (p(φ)) over positive parameters; for the oxygen extractionfraction E0 we use a uniform distribution over the interval[0, 1].These samples (φn, n = 1..N ) can then be used as anempirical approximation of the target distribution, e.g. forprediction of a BOLD signal y:p(y|D) =∫p(y|φ)p(φ|D)dφ≈ 1NN∑n=1p(y|φn), φn ∼ p(φ|D)The resulting approximate distribution clearly depends on theparticular data set D. To obtain information on the uncer-tainty of the posterior, we employ epoch-wise leave-one-outcross-validation. WithK quasi-independent epochs available,each split of the data leaves out one epoch for a test set, whichcan be used to validate the ability of the model to predictiontest data, while K − 1 are used for a training set to give anestimate of the posterior distribution.To get an estimate of the uncertainty of the posterior dis-tribution approximations, we fit a normal distribution to thesamples for each cross-validation split. This yields a distri-bution of means and variances that can be used to illustratethe variance of the distribution. Each split yields an unbiasedestimate of the true mean and variance of the distribution, andthe mean over all splits is a convex combination thus also anunbiased estimate, but with a reduced variance.953Authorized licensed use limited to: D anmarks Tekniske Informationscenter. D ownloaded on N ovember 28 , 2009 at 16:10  from IE E E  X plore .  R estrictions apply. 4. EXPERIMENTAL EVALUATIONThe described method was evaluated on a synthetically gen-erated data set and on a real BOLD fMRI data set.4.1. Synthetic dataThe synthetic data was obtained by simulating the hemody-namic model with parameters set to the maximum likelihoodvalues reported by Friston et al. [2],θ =[0.330.540.981.542.460.34 σ2w]T ,with σ2w set to produce a desired SNR (signal-to-noise ratio,measured as the ratio of the de-noised BOLD and observationnoise signals) close to 5.0 dB, which is similar to real record-ing conditions. κ and τi were set to 2.0 and 1.6 respectively.The model is initialized in x0, the states are evolved using aRunge-Kutta solver, and observations are made, adding Gaus-sian white noise with the prescribed variance. Each epochcontains 100 samples with sampling interval TR = 1.0s.To justify our assumption that the BOLD signal is inde-pendent between epochs, the stimulus for each epoch is set tozero for the last 30 seconds. This is helpful for two furtherreasons. Preprocessing (e.g. removing low-frequency noise),is aided in that such artifacts can be more accurately estimatedusing these ’resting’ portions of data. Finally, it allows us toassume a known, resting, physiological state (x0) at the startof each epoch. The stimulus is designed to evoke non-linearbehavior in the model; this is achieved by inter-stimulus in-tervals (ISI) and stimulus durations (SD) being sampled froma suitable gamma distribution.Figure 1 shows parts of the reconstructed neural and BOLDsignals. The posterior mean of κ and τ was 1.92 ± 0.51 and1.32 ± 0.22 respectively, resulting in very close reconstruc-tions.5 10 15 20 2500.10.20.30.40.50.60.70.80.91time(sec)neural activityTruePM(a)0 50 100 150 200 250!4!2024681012x 10!3time(sec)BOLD levelPMTrue(b)Fig. 1. (a): One pulse of the reconstructed and true neuralactivity. The shape of the signal is retrieved, although theslope is slightly off, corresponding to a bias in the τi estimate.(b): Reconstructed and true BOLD signal (first 3 epochs).4.2. Real dataThe data was acquired at Hvidovre Hospital, Denmark, us-ing a 3T scanner (Magnetom Trio, Siemens). We obtained1382 GRE EPI volumes each consisting of twelve 3mm slicesoriented along the calcarine sulcus. Additional parameterswhere TR=725 ms, TE=30 ms FOV=192 mm, 64x64 acqui-sition matrix, FA = 82. The stimulus consisted in a circularblack/white flickering checkerboard (24 degrees horizontal,18 degrees vertical) on a grey background. The checkers re-versed black/white at 8 Hz. The activation pattern (a(t)) usedto determine on- and offset of this stimulus was the same aswas used to generate the synthetic data.Fifty significantly activated (as determined by SPM2 anal-ysis 2 [8]) voxels in visual cortex were selected, and the meanof these was used as the BOLD signal.The results are shown in figure 2. The resulting shape ofthe neural model is similar to the one found for the syntheticdata, and is significantly different from square. The BOLD re-construction (see (3.1)) on test data is satisfactory, validatingthe model and method. We found posterior mean values of3.11± 0.76 for κ and 0.87± 0.19 (± one standard deviation)for τi; again, not close to a square pulse.Figure 3 shows the normal approximations to the posteriorhistograms. These indicate that although there is significantvariation, the mean of all the posterior means obtained fromthe leave-one-out cross-validation is identifiable. We expectwith longer sampling runs to bring down this variability - theimportant point is that this is a tool that gives guidance on thereliability of the estimated distribution.8 9 10 11 12 13 1400.20.40.60.81time(sec)neural activityPM+/! sd(a)0 50 100 150 200 250 300!505101520x 10!3time(sec)BOLD levelObservedPM(b)Fig. 2. (a): One pulse of the estimated neural activity. Theposterior mean reconstruction is shown together with the re-construction corresponding to± one standard deviation of thedistribution of the posterior mean. (b): Prediction on real testdata epochs (first 3 epochs).5. CONCLUSIONThe method detailed here can be used obtain to obtain approx-imate posterior distributions of model parameters together withestimates on the reliability of these approximations. For the2Software available from http://www.fil.ion.ucl.ac.uk/spm/954Authorized licensed use limited to: D anmarks Tekniske Informationscenter. D ownloaded on N ovember 28 , 2009 at 16:10  from IE E E  X plore .  R estrictions apply. 1 2 3 4 500.511.522.53κp(κ|D)(a)1 2 3 4 50123456κp(κ|D)(b)0 0.5 1 1.50123456τip(τ i|D)(c)0 0.5 1 1.5024681012τip(τ i|D)(d)Fig. 3. (a): Normal approximations of posterior parameterhistograms. (a) and (b) show the variation of the mean andvariance respectively, for the posterior distribution of κ. (c)and (d) show the same for τi (real data).parameters of the non-linear neural activity model, we foundposterior mean values of 3.11 ± 0.76 for κ and 0.87 ± 0.19(± one standard deviation) for τi. Both of these are on theedge of the ranges given in [3], although not statistically sig-nificantly so. There is little other information available onstatistical estimation of these parameters.We find that both κ and τ are identifiable for real BOLDdata, and that κ/τi for our data is significantly greater thanzero. Thus the square model of neural activity which is widelyused in BOLD analysis is not supported by our findings.The present model may be inverted to produce estimatesof neural activity as indicated in the work of Riera et al. [4]. In[4] a regularized radial basis function set is used, with param-eters estimated using a likelihood based approach which leadsto rather smooth activation estimates. Using our Bayesiansampling approach from an augmented posterior distributionincluding parameters of the neural activity time course (suchas stimulus onset times etc.) may be a way to let data de-termine the level of regularization, hence, potentially leadto more crisp estimates of non-trivial neural activation se-quences. This would be of particular interest in more complexactivation designs involving different stimulus activation con-ditions within epochs. In the present model we have focusedon the local hemodynamics in average data from a region.The BOLD hemodynamics is non-local and it is an impor-tant future task to produce a spatio-temporal hemodynamicmodel, which could also lead to improved spatial resolution.6. REFERENCES[1] R. B. Buxton, E. C. Wong, and L. R. Frank, “Dynamicsof blood flow and oxygenation changes during brain acti-vation: the balloon model,” MRM, vol. 39, pp. 855–864,June 1998.[2] K. J. Friston, A. Mechelli, R. Turner, and C.J. Price,“Nonlinear responses in fmri: the balloon model, volterrakernels, and other hemodynamics,” NeuroImage, vol. 12,pp. 466–477, 2000.[3] R. B. Buxton, K. Uludag, D. J. Dubowitz, and T. T. Liu,“Modeling the hemodynamic response to brain activa-tion,” NeuroImage, vol. 23, pp. 220–33, 2004.[4] J.J. Riera, J. Watanabe, K. Iwata, N. Miura, E. Aubert,T. Ozaki, and R. Kawashima, “A state-space model ofthe hemodynamic approach: nonlinear filtering of boldsignals,” Neuroimage, vol. 21, pp. 547–567, 2004.[5] T. Obata, T.T. Liu, K.L. Miller, W.M Luh, E.C. Wong,L.R. Frank, and R.B. Buxton, “Discrepancies betweenbold and flow dynamics in primary and supplementarymotor areas: application of the balloon model to the in-terpretation of bold transients,” Neuroimage, vol. 21, pp.144–153, 2004.[6] M. Abramowitz and I. A. Stegun, Eds., Handbookof Mathematical Functions with Formulas, Graphs, andMathematical Tables, 9th printing, New York: Dover,1972.[7] D.J.C. MacKay, Information theory, inference, and learn-ing algorithms, Cambridge University Press, 2003.[8] K.J. Friston, Introduction: Experimental design and sta-tistical parametric mapping, Academic Press, secondedition, 2003.955Authorized licensed use limited to: D anmarks Tekniske Informationscenter. D ownloaded on N ovember 28 , 2009 at 16:10  from IE E E  X plore .  R estrictions apply. 

Identification of non-linear models of neural activity in bold fmri

http://orbit.dtu.dk/ws/files/3964271/Jacobsen.pdf

Identification of non-linear models of neural activity in bold fmri

Abstract

Similar works

Full text

Available Versions

Crossref

Online Research Database In Technology