Abstract
Conventionally, as a preprocessing step, functional MRI (fMRI) data are spatially smoothed before further analysis, be it for activation mapping on task-based fMRI or functional connectivity analysis on resting-state fMRI data. When images are smoothed volumetrically, however, isotropic Gaussian kernels are generally used, which do not adapt to the underlying brain structure. Alternatively, cortical surface smoothing procedures provide the benefit of adapting the smoothing process to the underlying morphology, but require projecting volumetric data on to the surface. In this paper, leveraging principles from graph signal processing, we propose a volumetric spatial smoothing method that takes advantage of the gray-white and pial cortical surfaces, and as such, adapts the filtering process to the underlying morphological details at each point in the cortex.
I. INTRODUCTION
Functional magnetic resonance imaging (fMRI) is a key non-invasive imaging modality for studying brain activity based on the blood-oxygen-level-dependent signal [1]. When analyzing fMRI data, a conventional preprocessing step is to spatially smooth the data, which is generally done for two reasons: (i) to deal with the multiple testings through resorting to Gaussian random field theory [2]—as opposed to using stricter methods such as the Bonferroni correction—and (ii) to improve the data signal-to-noise ratio (SNR) by employing spatial smoothing such as lowpass filtering. With regards to (i), alternative methods have been proposed that obviate the need for Gaussian smoothing, namely, but not limited to, false discovery rate correction [3], permutation testing [4], and wavelet denoising [5]. With regards to (ii), the conventional use of isotropic filters—such as the Gaussian filter—comes at the expense of a loss in fine spatial details of the underlying activity [6], since in virtue of the matched filter argument, spatial filters are optimal in denoising data only if their shape and size conform to the underlying signal of interest. Given that the spatial profile of brain fMRI data is confined by the underlying neuroanatomical structure, alternative spatial filtering methods that incorporate knowledge about the underlying brain structure have been proposed, broadly classified into surface-based—see e.g. [7], [8], [9], [10], [11]—and volume-based—see e.g. [12], [13], [14], [15]—methods, which differ from methods such as [16], [17], [18], [19] in that they leverage an independent contrast image—different from the data to be smoothed—that confines the spatial profile of the filters.
In this paper, we propose a scheme that enables volumetric spatial smoothing of cerebral cortical fMRI data in such way that filtering is informed by the morphological structure of the cerebral cortex provided by cortical surface extractions. At the heart of the method lies the representation of a cerebral hemisphere cortex—a.k.a. the cortical ribbon—as a graph [20], [21], wherein each voxel is represented as a vertex and edges are defined based on the Euclidean adjacency of voxels but constrained by topological information provided by reconstructed white-gray border and pial cortical surfaces [22]. The graphs are subject-specific and thus encode the unique morphology of each individual’s cortical structure. Principles from the recently emerged field of graph signal processing [23], [24] are then leveraged to perform spatial filtering of fMRI data in such way that the filter profiles are adapted to cortical morphology.
II. METHODS
A. Graph signal processing: fundamentals
An undirected, non-weighted graph, denoted as , consists of a set of N vertices and a set of edges (i.e., pairs (i, j) where ), and can be represented by an N × N symmetric matrix A, i.e., the graph adjacency matrix, with its elements defined as ai,j = 1 if and ai,j = 0 otherwise. The graph’s symmetric normalized Laplacian matrix, denoted as L, is given as , where D denotes the graph’s degree matrix, which is a diagonal matrix with elements given as , and I denotes the identity matrix. Since L is real, symmetric, diagonally dominant, and with non-negative diagonal entries, it is positive semi-definite, and therefore, it can be eigen-decomposed as
(1) |
where Λ is a diagonal matrix that stores the eigenvalues λ1, …, λN ≔ λmax and U is a matrix that stores the corresponding eigenvectors in its columns, U = [u1|u2| ⋯ |uN]; hereon, we refer to eigenvectors of U as eigenmodes, in line with the convention in the neuroimaging community. The eigenvalues of U define a spectrum for [25], a space that can be seen as an extension of the Euclidean Fourier domain.
Given a graph , data residing on the vertices of can be seen as a graph signal, denoted as . Using U, f can be transformed to the spectrum , denoted as and obtained as , which is a representation that is commonly referred to as the graph Fourier transform of f. This transformation satisfies Parseval’s energy conservation relation [26], i.e.,. .
Given a graph and a graph signal f, it is beneficial to be able to implement filtering operations on f as in conventional signal processing. In particular, it is of interest to define graph filters that entail specific spectral properties, for instance, defining lowpass and highpass filters as in conventional signal processing. A graph filter can be conveniently defined in the graph Laplacian spectral domain as , using which, filtering of f, denoted as fh, can be straightforwardly implemented in the spectral domain by modulating with h(·), and then transforming the modulated signal to the vertex domain, given as
(2) |
where h(Λ) denotes a diagonal matrix with its entries computed by sampling h(·) at the eigenvalues. However, a shortcoming of this approach to graph signal filtering is that it requires Λ and U, i.e., eigen-decomposition of L. This is computationally highly impractical for large graphs, and, in particular, infeasable for the graphs proposed in this work that are typically of sizes in the range 100 K to 200 K vertices. An alternative approach is to approximate h(·) as a polynomial, denoted as , and to then implement filtering within the vertex domain as [27]
(3) |
(4) |
where in the last equality we used the property . As such, graph signal filtering can be seamlessly implemented through the mere use of polynomial matrix operations on L. In particular, it is insightful to observe that each row k in , where , is the vertex representation of spectral kernel h(·) when instantiated at vertex k; these set of vectors provide the impulse response set for spectral kernel h(·), which, unlike that in conventional signal processing, are shift-variant.
B. Dataset
We used the Human Connectome Project (HCP) dataset [28] in this study; the acquisition of the HCP dataset was approved by the Washington University Institutional Review Board and informed consent was obtained from all subjects. In particular, we used a subset of the dataset, the 100 unrelated subjects (54% female, mean age = 29.11± 3.67, age range = 22–36). To construct cortical graphs, we used the minimally preprocessed structural data, which come at 0.7 mm isotropic resolution, and the associated extended preprocessed version of the data, which included FreeSurfer surface extraction. We also used the HCP functional task fMRI data, which come at 2 mm isotopic resolution, and consisted of seven functional tasks: Emotion, Gambling, Language, Motor, Relational, Social, and Working Memory. The spatial smoothing methodology proposed in this work heavily relies on accurate co-registration between the structural and functional data, which is already implemented in the HCP preprocessed data. A full description of the imaging parameters and prepocessing steps can be found in [29].
C. Cerebral hemisphere cortex (CHC) graphs
We design subject-specific, voxel-resolution graphs that encode the morphology of CHC using the method presented in [20], [21], which is based on leveraging the volumetric representation of a subject’s CHC as extracted by FreeSurfer1 [22] from a T1-weighted MRI image, the so called ribbon representation. Each vertex of the 3D ribbon is considered as a graph vertex. Graph edges are preliminarily defined based on the adjacency of voxels within the cubic lattice of size 3×3×3, i.e., 26 neighborhood connectivity. Consequently, spurious edges that are anatomically unjustifiable—e.g. edges that connect voxels that lie on opposite sides of narrow sulci—are pruned out by using pial and white surface representations extracted using FreeSurfer. For a more detailed description of the design, we refer the interested reader to [20], [21].
D. Spectral graph heat kernel filters
Given a CHC graph, we design spatial smoothing filters associated with a heat kernel profile defined on the graph spectrum, i.e., k(λ) = e−τλ, ∀λ ∈ [0, λmax], where τ is a free parameter that specifies the resulting filter’s spatial extent; see Fig. 1. For a given τ, we approximated k(λ) by a Chebyshev polynomial—approximating a minimax polynomial, minimizing an upper bound on the approximation error [27]—and applied filtering as in (4).
E. Construction of cortical activity phantoms
A semi-synthetic phantom dataset was constructed through synthesizing cortical activation patterns and corrupting them with white Gaussian additive noise. The anatomical scans of the first 10 subjects from the HCP dataset were used to generate activation patterns that have confounding spatial shapes that are constrained by the cortical morphology of each individual subject, using a method similar to that proposed in [12]. Specifically, for each subject i, 10 random seed points were selected from within the cortical ribbon, treated as distinct centers of functional activity. The points were then represented as an indicator vector , wherein xi[n] = 1 if n corresponds to a center point and xi[n] = 0 otherwise, i.e., ∥xi∥1 = 10. An activation pattern that diffuses from the center points along the subject’s cortex was generated as yi = zi/ maxn{zi[n]}, where , and Ai denotes the adjacency matrix of subject i, and p is parameter that determines the extent of the diffused pattern. In all the following analyses, p was set to 5. The resulting activation patterns are unique to each subject and exhibit non-binary distributions of value within the range [0, 1]. The vector-formed activity patterns were then placed back into 3D volumetric representations; Fig. 2 shows the resulting phantoms across the 10 subjects. The phantoms were then corrupted with additive white Gaussian noise; phantoms with four noise levels were generated—standard deviations σn = 2, 4, 8 and 16, resulting in contrast-to-noise ratios (CNR)—maxn{yi}/σn [30]—of 1/2, 1/4, 1/8 and 1/16, respectively. Ten realizations were generated for each noise level, resulting in a total of 100 signals for each CNR.
III. RESULTS AND DISCUSSION
Fig. 3 shows realizations of a representative set of graph spatial smoothing filters when localized at different parts along the cortical ribbon. In contrast to conventional Gaussian smoothing wherein the spatial shape of the filtering kernel is the same at each position across the image, in graph-based spatial smoothing, the shape of graph filters adapts to each specific position along the domain they are realized at, based on the encoding provided by the graph.
To validate the strength of the proposed filters in denoising fMRI data, receiver operating characteristic (ROC) analysis was performed on the phantom set, see Section II-E. Each noise-added phantom was spatially smoothed with Gaussian spatial smoothing (GAUSS) and graph-based spatial smoothing (GRASS) over a range of filter sizes, FWHM = 1, 2, … 14 and τ = 10, 15, …, 100, respectively. For GAUSS, two approaches were compared: unconstrained GAUSS (uGAUSS) and constrained GAUSS (cGAUSS). In uGAUSS, the data were smoothed without masking out regions outside the ribbon; that is, the filter has access to values not only within the cortical ribbon but also to those in adjacent white matter or CSF. In cGAUSS, the region outside the ribbon was masked out (set to zero), and consequently, normalized convolution was performed; that is, filtering was done in three steps, first, the masked data was smoothed, then, a mask representing the cortical ribbon (with values one inside the ribbon and zero outside) was smoothed, and consequently, the former was divided by the latter.
After filtering, smoothed volumes were thresholded at 100 uniformly-spaced consecutive levels spanning the minimum and maximum value in each filtered volume to generate ROC curves—wherein the clean phantoms were treated as ground truth, and in turn, the area under curve of each ROC curve (AUCROC) was computed. Fig. 4 shows the results. Between the two GAUSS approaches, cGAUSS out-performed uGAUSS as it uses additional information about the delineation of the cortical ribbon, preventing mixture of pure noise from surrounding white matter and CSF. Across CNRs, GRASS showed superior performance over GAUSS. Although both cGAUSS and GRASS exploit cortical surface information, GRASS outperforms cGAUSS since i) it leverages filters whose spatial profile adapts to the domain, whereas in cGAUSS, the filters are merely masked by the cortical ribbon, and ii) it prevents mixing of values at touching parts of the cortex that are only adjacent in Euclidean sense but not in geodesic sense, for instance touching banks of narrow sulci. Best performance for GAUSS was obtained for FWHMs within the range 6 to 10, and for GRASS it was obtained for τ within the range 30 to 50.
We then performed Dice analysis to quantify the extent of difference between detection maps resulting from GRASS and cGAUSS on the real data. Given the large number of available combinations of cGAUSS and GRASS for different filter sizes, and lack of a one-to-one mapping between a Gaussian kernel of a given FWHM and a graph filter of a specific τ, we selected a subset of filter sizes for cGAUSS and GRASS for the analysis performed on the real data, based on the AUCROC performances obtained on the simulated data. In particular, we selected five filter sizes that provided the overall best performance for cGAUSS (FWHM = 6, 7, 8, 9 and 10) and GRASS (τ = 30, 35, 40, 45 and 50); see Fig. 4. With this selection of filter sizes, we are implicitly assuming that underlying activations in the real data have spatial spreads that are of approximately the same size as those in the simulated data. For each subject, and each functional task, the fMRI data were smoothed using the 10 above-listed filter settings, and activation mapping was then performed for one of the experimental conditions of each task (Emotion: fear, Gambling: win, Language: math, Motor: left hand, Relational: match, Social: mental, Working Memory: body 0-back; for details about the tasks and each experimental condition, we refer the interested reader to [31]). Resulting t-maps were then thresholded at 5% false discovery rate. For a given subject, a given task, and a given τ, the Dice similarity between the associated GRASS detection map and the five detection maps for cGAUSS were computed as dτ,fwhm = 2|Mτ ∩ Mfwhm|/|Mτ| + |Mfwhm|, for FWHM = 6, …, 10, where Mfwhm and Mτ denote the set of voxels in the detected map for cGAUSS and GRASS, respectively, and | · | denotes set cardinality. The maximum Dice was then selected, i.e., max{dτ,fwhm}fwhm=6, ⋯, 10. This analysis was repeated across subjects, tasks, and the five τ values. If activation mapping on data smoothed with neither cGAUSS nor GRASS resulted in any detections, Dice similarity was not computed.
Fig. 5(a)-left shows the maximum Dice for different tasks for τ = 40, showing values within the range 0 to 0.92 and an average of 0.8 across the seven tasks2. A similar average Dice similarity was observed using the other four τ values; see Fig. 5(a)-right. Despite the similarity of detection maps for GRASS and cGAUSS, the two approaches nevertheless resulted in a notable set of detections that are not common between the two. Figs. 5(b) and (c) show the extent of unique detections by cGAUSS and GRASS, respectively, as ratios of total number of detections by each method; in particular, Fig. 5(b) shows |Mfwhm| − |Mτ ∩ Mfwhm|, and Fig. 5(c) shows |Mτ| − |Mτ ∩ Mfwhm|, for the τ and the FWHM for which Dice values are shown in Fig. 5(a). On average, 21%3 of the detections by cGAUSS were not present in the detection maps for GRASS, and 17%4 vice versa. Due to lack of ground truth, it is not possible to make any concrete statements when interpreting these results. However, if we are to accept the validity of the simulation results—specifically, that the generated synthetic phantoms resemble true activations—the superior performance of GRASS over cGAUSS on the simulations suggests the potential of having obtained a better specificity and sensitivity using GRASS.
IV. CONCLUSIONS
We proposed a method that enables volumetric spatial smoothing of functional data in such way that filtering is informed by the morphological structure of the cerebral cortex, as defined by extracted cortical surfaces. In future work, we will compare the proposed method with spatial smoothing performed on cortical surface [10], [14], and will also explore the benefits of the proposed method for processing high spatial resolution fMRI data [32], [33].
Acknowledgments
This work, and HB, were supported by the Swedish Research Council (2018-06689), and in part by the Royal Physiographic Society of Lund, the Thorsten and Elsa Segerfalk Foundation, the Hans Werthén Foundation and the Sweden-America Foundation. CFW was supported by the NIH grants P41EB015902 and R01MH074794. IA was supported by the NIH grants K01DK101631 and R56AG068261.
Footnotes
An implementation of the proposed method will be made available as a MATLAB package at https://github.com/aitchbi/GRASS.
Emotion: 0.8, Gambling: 0.84, Language: 0.74, Motor: 0.76, Relational: 0.81, Social: 0.85 and WM: 0.8.
Emotion: 22%, Gambling: 15%, Language: 26%, Motor: 29%, Relational: 18%, Social:14% and WM: 21%.
Emotion: 15%, Gambling: 14%, Language: 24%, Motor: 19%, Relational: 13%, Social: 16% and WM: 16%.
References
- [1].Ogawa S, Menon RS, Tank DW, Kim SG, Merkle H, Ellermann JM, and Ugurbil K, “Functional brain mapping by blood oxygenation level-dependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model.,” Biophys. J, vol. 64, no. 3, pp. 803–812, 1993. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [2].Poline JB, Worsley KJ, Evans AC, and Friston KJ, “Combining spatial extent and peak intensity to test for activations in functional imaging.,” Neuroimage, vol. 5, pp. 83–96, 1997. [DOI] [PubMed] [Google Scholar]
- [3].Genovese CR, Lazar NA, and Nichols T, “Thresholding of statistical maps in functional neuroimaging using the false discovery rate,” Neuroimage, vol. 15, no. 4, pp. 870–878, 2002. [DOI] [PubMed] [Google Scholar]
- [4].Eklund A, Andersson M, and Knutsson H, “Fast random permutation tests enable objective evaluation of methods for single-subject fMRI analysis,” Int. J. Biomed. Imag, vol. 2011, 2011. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [5].Van De Ville D, Blu T, and Unser M, “Integrated wavelet processing and spatial statistical testing of fMRI data,” Neuroimage, vol. 23, pp. 1472–1485, 2004. [DOI] [PubMed] [Google Scholar]
- [6].Geissler A, Lanzenberger R, Barth M, Tahamtan AR, Milakara D, Gartus A, and Beisteiner R, “Influence of fMRI smoothing procedures on replicability of fine scale motor localization,” Neuroimage, vol. 24, no. 2, pp. 323–331, 2005. [DOI] [PubMed] [Google Scholar]
- [7].Kiebel SJ, Goebel R, and Friston KJ, “Anatomically informed basis functions,” Neuroimage, vol. 11, no. 6, pp. 656–667, 2000. [DOI] [PubMed] [Google Scholar]
- [8].Chung MK and Taylor J, “Diffusion smoothing on brain surface via finite element method,” in Proc. IEEE Int. Symp. Biomed. Imaging IEEE, 2004, pp. 432–435. [Google Scholar]
- [9].Qiu A, Bitouk D, and Miller MI, “Smooth functional and structural maps on the neocortex via orthonormal bases of the Laplace-Beltrami operator,” IEEE Trans. Med. Imag, vol. 25, no. 10, pp. 1296–1306, 2006. [DOI] [PubMed] [Google Scholar]
- [10].Blazejewska AI, Fischl B, Wald LL, and Polimeni JR, “Intracortical smoothing of small-voxel fMRI data can provide increased detection power without spatial resolution losses compared to conventional large-voxel fMRI data,” Neuroimage, vol. 189, pp. 601–614, 2019. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [11].Huang SG, Lyu I, Qiu A, and Chung MK, “Fast polynomial approximation of heat kernel convolution on manifolds and its application to brain sulcal and gyral graph pattern analysis,” IEEE Trans. Med. Imag, vol. 39, no. 6, pp. 2201–2212, 2020. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [12].Behjat H, Leonardi N, Sörnmo L, and Van De Ville D, “Anatomically-adapted graph wavelets for improved group-level fMRI activation mapping,” Neuroimage, vol. 123, pp. 185–199, 2015. [DOI] [PubMed] [Google Scholar]
- [13].Behjat Hamid and Van De Ville Dimitri, “Spectral design of signal-adapted tight frames on graphs,” in Vertex-Frequency Analysis of Graph Signals, pp. 177–206. Springer, 2019. [Google Scholar]
- [14].Huber LR, Poser BA, Bandettini PA, Arora K, Wagstyl K, Cho S, Goense J, Nothnagel N, Morgan AT, Mueller AK, et al. , “LayNii: A software suite for layer-fMRI. revision 1,” BioRxiv, 2020. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [15].Abramian D, Larsson M, Eklund A, Aganj I, Westin CF, and Behjat H, “Diffusion-informed spatial smoothing of fMRI data in white matter using spectral graph filters,” Neuroimage, in press, 2021. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [16].Smith SM and Brady JM, “SUSAN — a new approach to low level image processing,” Int. J. Comput. Vis, vol. 23, no. 1, pp. 45–78, 1997. [Google Scholar]
- [17].Friman O, Borga M, Lundberg P, and Knutsson H, “Adaptive analysis of fMRI data,” Neuroimage, vol. 19, no. 3, pp. 837–845, 2003. [DOI] [PubMed] [Google Scholar]
- [18].Rydell J, Knutsson H, and Borga M, “Bilateral filtering of fMRI data,” IEEE J. Sel. Top. Signal Process, vol. 2, no. 6, pp. 891–896, 2008. [Google Scholar]
- [19].Lohmann G, Stelzer J, Lacosse E, Kumar VJ, Mueller K, Kuehn E, Grodd W, and Scheffler K, “LISA improves statistical analysis for fMRI,” Nat. Commun, vol. 9, no. 1, pp. 1–9, 2018. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [20].Behjat H and Larsson M, “Spectral characterization of functional MRI data on voxel-resolution cortical graphs,” in Proc. IEEE Int. Symp. Biomed. Imaging, 2020, pp. 558–562. [Google Scholar]
- [21].Maghsadhagh S, Eklund A, and Behjat H, “Graph spectral characterization of brain cortical morphology,” in Proc. IEEE Int. Conf. Eng. Med. Biol. Soc, July 2019, pp. 458–462. [DOI] [PubMed] [Google Scholar]
- [22].Fischl B, “Freesurfer,” Neuroimage, vol. 62, no. 2, pp. 774–781, 2012. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [23].Shuman DI, Narang SK, Frossard P, Ortega A, and Vandergheynst P, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” IEEE Signal Process. Mag, vol. 30, no. 3, pp. 83–98, 2013. [Google Scholar]
- [24].Ortega A, Frossard P, Kovačević J, Moura JMF, and Vandergheynst P, “Graph signal processing: Overview, challenges, and applications,” Proc. IEEE, vol. 106, no. 5, pp. 808–828, May 2018. [Google Scholar]
- [25].Chung FRK, Spectral graph theory, AMS, Providence, RI, 1997. [Google Scholar]
- [26].Shuman DI, Ricaud B, and Vandergheynst P, “Vertex-frequency analysis on graphs,” Appl. Comput. Harmon. Anal, 2015, doi: 10.1016/j.acha.2015.02.005. [DOI] [Google Scholar]
- [27].Hammond DK, Vandergheynst P, and Gribonval R, “Wavelets on graphs via spectral graph theory,” Appl. Comput. Harmon. Anal, vol. 30, no. 2, pp. 129–150, 2011. [Google Scholar]
- [28].Van Essen DC, Smith SM, Barch DM, Behrens TE, Yacoub E, Ugurbil K, and WU-Minn HCP Consortium., “The WU-Minn human connectome project: an overview.,” Neuroimage, vol. 80, pp. 62–79, 2013. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [29].Glasser Matthew F, Sotiropoulos Stamatios N, Wilson J Anthony, Coalson Timothy S, Fischl Bruce, Andersson Jesper L, Xu Junqian, Jbabdi Saad, Webster Matthew, Polimeni Jonathan R, et al. , “The minimal preprocessing pipelines for the human connectome project,” Neuroimage, vol. 80, pp. 105–124, 2013. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [30].Welvaert M and Rosseel Y, “On the definition of signal-to-noise ratio and contrast-to-noise ratio for fMRI data,” PloS one, vol. 8, no. 11, pp. e77089, 2013. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [31].HCP WU-Minn, “1200 subjects data release reference manual,” URL https://www.humanconnectome.org, 2017.
- [32].Moerel M, De Martino F, Uğurbil K, Yacoub E, and Formisano E, “Processing complexity increases in superficial layers of human primary auditory cortex,” Sci. Rep, vol. 9, no. 1, pp. 1–9, 2019. [DOI] [PMC free article] [PubMed] [Google Scholar]
- [33].Finn ES, Huber L, Jangraw DC, Molfese PJ, and Bandettini PA, “Layer-dependent activity in human prefrontal cortex during working memory,” Nat. Neurosci, vol. 22, no. 10, pp. 1687–1695, 2019. [DOI] [PMC free article] [PubMed] [Google Scholar]