We propose a new method for the detection of activated voxels in event-related BOLD fMRI data. We model the statistics of the wavelet histograms derived from each voxel time series independently through a generalized Gaussian distribution (GGD). We perform k-means clustering of the GGDs characterizing the voxel data in a synthetic data set, using the symmetrized Kullback-Leibler divergence (KLD) as a similarity measure. We compare our technique with GLM modeling and with another clustering method for activation detection that directly uses the wavelet coefficients as features. Our method is shown to be considerably more stable against realistic hemodynamic variability

Rosseel, Yves

Verdoolaege, Geert

Ghent University Academic Bibliography

ACTIVATION DETECTION IN EVENT-RELATED FMRITHROUGH CLUSTERING OF WAVELET DISTRIBUTIONSGeert Verdoolaege and Yves RosseelDepartment of Data Analysis, Ghent University,Henri Dunantlaan 1, 9000 Gent, BelgiumABSTRACTWe propose a new method for the detection of activated vox-els in event-related BOLD fMRI data. We model the statis-tics of the wavelet histograms derived from each voxel timeseries independently through a generalized Gaussian distri-bution (GGD). We perform k-means clustering of the GGDscharacterizing the voxel data in a synthetic data set, using thesymmetrized Kullback-Leibler divergence (KLD) as a simi-larity measure. We compare our technique with GLM mod-eling and with another clustering method for activation de-tection that directly uses the wavelet coefficients as features.Our method is shown to be considerably more stable againstrealistic hemodynamic variability.Index Terms— fMRI, k-means clustering, generalizedGaussian distribution, Kullback-Leibler divergence1. INTRODUCTIONIn neuroscience, functional magnetic resonance imaging(fMRI) has become one of the most widely used techniquesfor the mapping of brain activity that is related to some cog-nitive task. However, various sources of noise are presentin the measured BOLD signal and the signal-to-noise ratio(SNR) in activated areas is typically low. This renders thedetection of activated voxels a challenging task and manydata analysis techniques, both hypothesis- and data-driven,have been proposed in the past to accomplish this.The BOLD response is often modeled through the convo-lution of a neural activation signal with a hemodynamic re-sponse function (HRF). The general linear model (GLM) is awell-known tool in statistics that is frequently used for fMRIactivation detection [1]. The GLM basically compares theacquired BOLD time series to a predefined set of basis func-tions. In many cases it can be useful to incorporate some apriori knowledge about the expected hemodynamic responseinto the basis functions, hence constructing a so-called in-formed basis. However, it has been shown that there can bea considerable variability in the HRF from one brain regionto another, or even in the same region on a trial-by-trial ba-sis [2]. In such a case it may be better to use a more genericset of basis functions and in this respect the wavelet basis isparticularly well suited. Indeed, it has been demonstrated thatwavelets provide a natural framework for the analysis of neu-robiological signals, which often exhibit fractal properties andnonstationary behavior of interest at multiple scales (see [3],which discusses other application of wavelets in fMRI dataanalysis as well). This provides our motivation for character-izing fMRI data through features in the wavelet domain.Although many fMRI data analysis schemes, such as theGLM, operate on a voxelwise basis, there are two main disad-vantages attached to this mode of operation. First, conductingstatistical tests in a mass-univariate way poses a nontrivialmultiple comparisons problem. Second, a voxelwise methoddoes not take into account any spatial correlation structurewithin the brain and is unable to detect activation patterns.In this paper, we use a data-driven clustering approach ofevent-related fMRI time series, which does not suffer fromthe above drawbacks. Cluster analysis of fMRI data hasbeen considered before, but the employed features are oftendirectly linked to the raw time series or they depend on aspecific HRF model (see [4] and references therein). Hence,such features are in general less robust against HRF variabil-ity, as we will show in the present work. For instance, in [5]clustering was performed on a subset of the detail coefficientsresulting from a discrete wavelet transform of the voxel timeseries (which we will refer to as the method ‘COEF’). Theclustering method that we propose, which we here call themethod ‘STAT’, is conceptually simple and is also basedon wavelet features. However, instead of clustering directlyon the wavelet coefficients, we use the statistics of waveletdetail coefficients as an indicator for activation. We modelthe wavelet statistics via a generalized Gaussian distribution(GGD) and we evaluate GGD similarity through the sym-metrized Kullback-Leibler divergence (KLD). The KLD issubsequently applied in a k-means clustering algorithm. Ina comparison of our technique with both COEF and classicGLM modeling, we show that the use of wavelet statisticsleads to an enhanced stability of the method against hemody-namic variability. We work with an artificial data set, whichhas the advantage that the activation pattern is known a pri-ori, facilitating comparison of different activation detectionstrategies.This contribution is organized as follows. In Section 2 wepresent the design of our experiments: the generation of theartificial data set, the statistics of the wavelet coefficients andGGD similarity measurement via the KLD in a k-means clus-tering scheme. In Section 3 we discuss the results of the clus-tering experiments for activation detection using STAT andwe compare with the GLM and COEF. Section 4 provides aconclusion and an outlook towards future work.2. DESIGN OF CLUSTERING EXPERIMENTS2.1. Synthetic data setWe demonstrate activation detection on an artificial data setconsisting of a volume of 14 × 14 × 14 voxels, every voxelcontaining a time series of 336 samples. A repetition time(TR) of 2 s was chosen. Activation was introduced in half ofthe number of time series (voxels), consisting of 24 regularlytimed block stimuli of 4 s duration each. The BOLD responsecorresponding to this activation was obtained by convolvingwith the HRF. In a first experiment a fixed canonical HRFthroughout the simulated volume was assumed. The canoni-cal HRF h(t) provides a physiologically reasonable model fora typical hemodynamic response [2]. It is defined as follows:h(t) =(td1)a1exp(d1 − tb1)−c(td2)a2exp(d2 − tb2).Here, d1 = a1b1 is the time associated to the maximum ofthe initial peak and d2 = a2b2 is the time to the minimumof the undershoot. We started from typical values for theparameters: a1 = 6, a2 = 12, b ≡ b1 = b2 = 0.9 andc = 0.35. The corresponding function h(t) is drawn in Fig-ure 1. A second experiment was designed in order to studythe robustness of the activation detection scheme against spa-tial hemodynamic variability. To this end a fluctuation wasintroduced in the HRF dispersion (parameter b) and ampli-tude of the undershoot (c). In addition, a variable responselatency was imposed by delaying the HRF by a time lag D.The values of b, c and D were chosen at random from voxelto voxel, by sampling from the following truncated Gaussiandistributions:b : N (0.9, 1.0) ∧ 0.5 ≤ b ≤ 2.0,c : N (0.35, 1.0) ∧ 0.0 ≤ c ≤ 1.0,D : N (0.0, 1.5) ∧ 0.0 ≤ D ≤ 2.0,where N (µ, σ) signifies the normal distribution with mean µand standard deviation σ. This represents a realistic hemody-namic variability supported by empirical data [2]. Apart fromthe canonical HRF, some more example HRFs are drawn inFigure 1. No trial-to-trial HRF variability was assumed.Several models of varying complexity for the noise infMRI data are commonly used. In the present experiments0 5 10 15 20 25 30 35 40−0.8−0.400.40.8Time (s)HRF (A.U.)  b = 0.9, c = 0.35 (canonical)b = 0.5, c = 0.35b = 1.5, c = 0.35b = 0.9, c = 0.10b = 0.9, c = 0.80Fig. 1. Canonical HRF and examples of other HRFs for someextreme values of the parameters in the HRF model that werevaried from voxel to voxel.we simply assumed additive Gaussian noise leading to a rel-atively low SNR. The (functional) SNR in fMRI is usuallymentioned as the ratio of the standard deviation σs of the puresignal (activation), to the standard deviation σn of the noise.Typical values are σs/σn . 1, while we assumed an SNRof 0.6 (approximately −4 dB). The non-activated time seriescontained only noise with the same standard deviation. Anexample of a simulated time series, showing both the pureand the noisy signal, is given in Figure 2.2.2. Wavelet statisticsApplied to our synthetic data set, the purpose of an activa-tion detection scheme based on clustering is to discriminatebetween the activated half of the simulated volume and thenon-activated part. The selection of features used for ourclustering method STAT was conducted as follows. First, theredundant discrete wavelet transform (RDWT) was appliedindependently to each time series, using Daubechies’ 4-tapwavelets. An example of the RDWT of an artificial activatedtime series (BOLD response + noise) is presented in Fig-ure 2, showing the first four levels. The effect of the originalBOLD response can easily be discerned at the higher wave-let detail scales. The RDWT results in a number of wave-let detail coefficients sufficient to adequately characterize thewavelet statistics. Wavelet detail histograms, which have zeromean, often exhibit a distinctly non-Gaussian shape. For in-stance, in [6] wavelet coefficients for gray-level images weremodeled in an image retrieval context using a univariate zero-mean generalized Gaussian distribution, also known by thename of exponential power distribution. The GGD in generalprovided a better fit than the Gaussian density. The univari-ate zero-mean GGD is defined through the following densityfunction:p(x|α, β) =β2αΓ(1/β)exp[−(|x|α)β],00.51BOLD(A.U.)−202BOLD+ noise(A.U.)−202Detailscale 1−202Detailscale 2−202Detailscale 30 100 200 300 400 500 600 700−202Time (s)Detailscale 4Fig. 2. The BOLD response (top frame) as a result of a syn-thetic regular series of block stimuli (4 s each), employingthe canonical HRF. The second frame shows the same signalwith added Gaussian noise (−4 dB). In the following framesthe wavelet detail coefficients are plotted corresponding to thefirst four levels of the RDWT of the noisy signal using theDaubechies 4-tap wavelets.where Γ(.) represents the Gamma function. α models thewidth of the density function (related to the standard devia-tion), while β is connected to its kurtosis. Clearly, the Gaus-sian distribution is obtained by setting β = 2.Example histograms of the redundant wavelet detail coef-ficients at the fourth scale obtained from both a non-activatedand an activated times series (same data as in Figure 2), aredisplayed in Figure 3. The respective GGD fits as well as thebest-fit Gaussian distributions are also shown and it can beseen that the GGD captures the wavelet statistics much bet-ter than the Gaussian. Accordingly, we described the waveletdetail coefficients corresponding to the time series in our ar-tificial volume at the first four levels independently through aset of univariate GGDs. The choice of wavelet scales largelydepends on the nature of the signal, but in our applicationwe found that including the first scale as well as larger scalesthan the fourth had little impact on clustering performance,so we ignored them in our detection scheme. This measureis, in addition, beneficial for the COEF method since it sub-stantially limits the number of features. The features that weused to characterize each voxel time series are thus given bythe parameters α and β of the univariate GGDs modeling thewavelet statistics. Since we routinely employed three waveletscales (scales two to four) we obtained six features per voxel,representing a dramatic reduction of dimensionality.2.3. ClusteringThe next step in our STAT method involved the clustering ofvoxels on the basis of their corresponding GGD features. Anoverview and comparison of several clustering algorithms inthe context of fMRI cluster analysis was given in [4]. The su-−4 −3 −2 −1 0 1 2 3 400.20.40.60.81Wavelet coefficientsRelative frequency  Non−activatedActivatedFig. 3. Histograms for the redundant wavelet detail coeffi-cients at level four corresponding to a time series containingonly noise (non-activated) and a time series in an activatedvoxel. The GGD fits are displayed as well (full lines) togetherwith the best-fit Gaussian density functions (dashed lines).perior clustering performance of the k-means and neural gasalgorithms was demonstrated, compared to other well-knowntechniques. For this reason in the present study we appliedthe k-means algorithm, which in addition is simple to imple-ment. The k-means algorithm assumes the number of clustersk is known and represents the clusters by their centroid. Theaim is to minimize the within-cluster sum-of-squares. In ourcase k = 2 since we simply distinguish between activatedand non-activated voxels. The initialization phase of the al-gorithm consists of assigning random cluster centroids. Thealgorithm next proceeds iteratively by cycling through twosteps: first, assign each point to the cluster with the nearestcentroid. Then, for every obtained cluster recalculate the cen-troid. The clustering efficiency depends partly on the choiceof similarity measure, for which we took the symmetrizedKullback-Leibler divergence. A closed-form expression ex-ists for the KLD between univariate zero-mean GGDs [6],rendering our implementation of the k-means algorithm suf-ficiently fast. In addition, the wavelet scales were assumed tobe independent, so that the total KLD between voxels couldsimply be obtained by summing the KLDs between corre-sponding scales [6]. The cluster centroids were defined bythe Karcher mean of the cluster points, again using the sym-metrized KLD as a similarity measure.3. RESULTS AND DISCUSSIONAs mentioned in the introduction, we compared our techniqueto the method COEF, for which we used as features the wave-let detail coefficients in scales two to four obtained with theDWT (including subsampling). In COEF the Euclidean dis-tance between wavelet coefficients is applied in a k-meansclustering algorithm. In both the case of COEF and STAT, weperformed the k-means cluster analysis several times in or-der to study the effect of random initialization of cluster cen-troids, which turned out to be negligible. A typical set of clus-ters in GGD parameter space is shown in Figure 4, together1 1.2 1.4 1.6123αβ  1 1.2 1.4 1.6 1.81234αβ1 1.4 1.8 2.212345αβ1 1.4 1.8 2.202468αβ  Non−activated voxelNon−activated cluster centerActivated voxelActivated cluster centerScale 1 Scale 2Scale 4Scale 3Fig. 4. GGD parameters, displayed independently for the firstfour scales of the RDWT, corresponding to the voxel timeseries in the synthetic volume (each point represents a voxel).Cluster centroids obtained via a k-means clustering algorithmare also indicated.with the cluster centroids found by the k-means algorithm. Aclear non-Gaussian behavior can be observed of the waveletstatistics in many voxels. It should be noted that the clustershapes do not necessarily reflect the true similarity structureof the voxels, since the plot is in Euclidean space while theapplied similarity measure is non-Euclidean (KLD).The results from the two experiments (fixed and variableHRF) are summarized in Table 1 in terms of false positiverates (FPR) and true positive rates (TPR). For comparisonwith the performance of the GLM, using the canonical HRFas a model, the TPR for a given FPR was determined from thedistributions of the test statistic for the non-activated and ac-tivated voxels. Since for both STAT and COEF the FPR turnsout to lie below 1%, we chose to mention the TPR obtainedby the GLM corresponding to an FPR of 0.1% and 1%.All three methods perform very similarly in case of a fixedcanonical HRF throughout the simulated volume. However,when HRF variability is introduced, the FPRs remain stablebut our STAT method exhibits a substantially higher TPR thanboth the GLM and the COEF method. Hence, STAT turns outto be much more robust against HRF variability compared tothe GLM and COEF.4. CONCLUSIONWe have proposed a new data-driven method for activationdetection in event-related BOLD fMRI. The technique isbased on k-means clustering of the parameters of a general-ized Gaussian distribution, fit to wavelet detail coefficientsobtained from the fMRI time series via the redundant discretewavelet transform. We have compared our method to classicGLM modeling and to a clustering-based scheme that usesraw wavelet coefficients as features. Our method turns out toGLM COEF STATCanonical HRF FPR (%) 0.1 / 1.0 0.1 0.3TPR (%) 100 100 100Variable HRF FPR (%) 0.1 / 1.0 0.0 0.0TPR (%) 83.3 / 88.1 54.7 92.5Table 1. False positive rates (FPR) and true positive rates(TPR) in the case of a fixed and variable HRF in an experi-ment on an artificial fMRI volume using the methods GLM,COEF and STAT.be considerably more stable against hemodynamic variability,which is ubiquitous in real fMRI data. A similar robustnessagainst uncertainty in stimulus timing is to be expected (use-ful for e.g. resting-state experiments). This will be the subjectof further study, as well as the use of other clustering algo-rithms (e.g. Gaussian mixture modeling) and GGD similaritymeasures (e.g. the geodesic distance).5. REFERENCES[1] K.J. Friston, A.P. Holmes, K.J. Worsley, J.-B. Poline,C.D. Frith, and R.S.J. Frackowiak, “Statistical ParametricMaps in functional imaging: A general linear approach,”Hum. Brain Mapp., vol. 2, pp. 189–210, 1994.[2] K.J. Friston, P. Fletcher, O. Josephs, A. Holmes, M.D.Rugg, and R. Turner, “Event-related fMRI: Characteriz-ing differential responses,” Neuroimage, vol. 7, pp. 30–40, 1998.[3] E. Bullmore, J. Fadili, V. Maxim, L. Xendur, B. Whitcher,J. Suckling, M. Brammer, and M. Breakspear, “Waveletsand functional magnetic resonance imaging of the humanbrain,” Neuroimage, vol. 23, pp. S234–S249, 2004.[4] E. Dimitriadou, M. Barth, C. Windischberger, K. Hornik,and E. Moser, “A quantitative comparison of functionalMRI cluster analysis,” Artif. Intell. Med., vol. 31, pp. 57–71, 2004.[5] B. Whitcher, A.J. Schwarz, H. Barjat, S.C. Smart, R.I.Grundy, and M.F. James, “Wavelet-based cluster anal-ysis: data-driven grouping of voxel time courses withapplication to perfusion-weighted and pharmacologicalMRI of the rat brain,” Neuroimage, vol. 24, pp. 281–295,2005.[6] M.N. Do and M. Vetterli, “Wavelet-based texture re-trieval using generalized Gaussian density and Kullback-Leibler distance,” IEEE Trans. Image Process., vol. 11,pp. 146–158, 2002.

Activation detection in event-related FMRI through clustering of wavelet distributions

Geert Verdoolaege

Yves Rosseel

Crossref

Activation detection in event-related fMRI through clustering ofwavelet distributions

https://biblio.ugent.be/publication/1234913/file/6752641.pdf

Activation detection in event-related FMRI through clustering of wavelet distributions

Abstract

Similar works

Full text

Available Versions

Ghent University Academic Bibliography

Crossref