# Improved functional MRI activation mapping in white matter through
diffusion-adapted spatial filtering

###### Abstract

Brain activation mapping using functional MRI (fMRI) based on blood
oxygenation level-dependent (BOLD) contrast has been conventionally focused on
probing gray matter, the BOLD contrast in white matter having been generally
disregarded. Recent results have provided evidence of the functional
significance of the white matter BOLD signal, showing at the same time that
its correlation structure is highly anisotropic, and related to the diffusion
tensor in shape and orientation. This evidence suggests that conventional
isotropic Gaussian filters are inadequate for denoising white matter fMRI
data, since they are incapable of adapting to the complex anisotropic domain
of white matter axonal connections. In this paper we explore a graph-based
description of the white matter developed from diffusion MRI data, which is
capable of encoding the anisotropy of the domain. Based on this representation
we design localized spatial filters that adapt to white matter structure by
leveraging graph signal processing principles. The performance of the proposed
filtering technique is evaluated on semi-synthetic data, where it shows
potential for greater sensitivity and specificity in white matter activation
mapping, compared to isotropic filtering.

Index Terms—  functional MRI, diffusion MRI, white matter, adaptive filtering

##  1 Introduction

Functional magnetic resonance imaging (fMRI) is a noninvasive technique that
allows investigation of the activity of the brain while performing a task or
at rest. Most fMRI studies report activations in the gray matter, while
reports of white matter activations are relatively sparse [1, 2, 3, 4, 5]. The
significance of such reports is controversial [6], as the sources of the
blood-oxygenation level dependent (BOLD) signal in white matter are not fully
understood. Despite this, significant evidence directly linking the BOLD
signal in white matter to neural activity has recently been presented [7].

The relative scarcity of reports of white matter activations can be partially
explained by the anatomical and physiological differences between white and
gray matter [6], which suggest the potential need for different experimental
designs and analysis methods to optimize detection power in each tissue type.
It has been shown, for example, that increased T2-weighting improves the
sensitivity to callosal activations in an interhemispheric transfer task [8].
Therefore, we expect the development of methods geared specifically towards
white matter to be required in order to investigate the functional
significance of the white matter BOLD signal.

An important distinguishing feature of the BOLD signal in white matter is that
it exhibits a spatial correlation structure grossly consistent with the
directions of water diffusion, as measured by diffusion tensor imaging (DTI)
[9]. This correlation structure is present during rest and becomes even more
pronounced under functional loading [10, 7]. As a consequence, isotropic
Gaussian filtering, commonly used to increase the signal-to-noise ratio of the
fMRI signal, may prove inadequate for use on the highly anisotropic white
matter domain.

The anisotropy of the white matter domain makes it suited for a graph-based
description. Graph signal processing methods have recently experienced
significant development [11, 12], resulting in their successful application to
the analysis of the BOLD signal in gray matter [13, 14]. In this work we
propose a graph-based adaptive filtering approach for the white matter BOLD
signal which incorporates information from diffusion MRI in order to adapt the
filter design to the white matter mircostructure. We test our proposed
filtering approach on semi-synthetic phantoms and show its increased
sensitivity and specificity for activations in white matter.

##  2 Methods

###  2.1 Preliminaries

We consider undirected, connected, weighted graphs
$\mathcal{G}=(\mathcal{V},\mathcal{E},\mathbf{A})$, where $\mathcal{V}$ is a
set of vertices with $|\mathcal{V}|=N_{v}$, $\mathcal{E}$ a set of edges
connecting pairs $(i,j)$ of vertices, and $\mathbf{A}$ is an adjacency matrix
whose nonzero elements $a_{i,j}$ represent the weight of edges
$(i,j)\in\mathcal{E}$. The diagonal degree matrix, denoted $\mathbf{D}$,
associated to $\mathbf{A}$ is defined with elements $d_{i,i}=\sum_{j}a_{i,j}$.

The normalized Laplacian matrix of $\mathcal{G}$ can then be obtained as
$\mathcal{L}=\mathbf{I}-\mathbf{D}^{-1/2}\mathbf{A}\mathbf{D}^{-1/2}$. As
$\mathcal{L}$ is a real symmetric matrix, it can be diagonalized, resulting in
a set of $N_{v}$ real non-negative eigenvalues which define the graph spectrum
$\mathbf{\Lambda}$, i.e.,
$\mathbf{\Lambda}=\\{0=\lambda_{1}\leq\lambda_{2}\ldots\leq\lambda_{N_{v}}\stackrel{{\scriptstyle\text{def}}}{{=}}\lambda_{\text{max}}\leq
2\\}$. The associated eigenvectors, denoted
$\\{\bm{\chi}_{l}\\}_{l=1}^{N_{v}}$, form an orthonormal basis spanning the
$\ell^{2}(\mathcal{G})$ space. Here $\ell^{2}(\mathcal{G})$ represents the
Hilbert space of all square-integrable graph signals
$\bm{f}:\mathcal{V}\to\mathbb{R}$ defined on the vertex set $\mathcal{V}$. A
graph signal can therefore be seen as a vector
$\bm{f}\in\ell^{2}(\mathcal{G})$ whose $n$-th component represent the signal
value at the $n$-th vertex of $\mathcal{G}$. A graph signal $\bm{f}$ can be
filtered with spectral kernel $\mathcal{K}(\lambda)$ as

$\left(F_{\bm{k}}\bm{f}\right)[m]=\sum_{l=1}^{N_{v}}\hat{\bm{k}}[l]\hat{\bm{f}}[l]\bm{\chi}_{l}[m],$
(1)

where $\hat{\bm{k}}$ is the sampled version of $\mathcal{K}(\lambda)$ obtained
as $\hat{\bm{k}}[l]=\mathcal{K}(\lambda_{l}),l=1,\dots,N_{v}$ and
$\hat{\bm{f}}$ is the graph Fourier transform of $\bm{f}$ (see [15] for an
introduction to graph signal processing).

###  2.2 White matter graph construction

We define a white matter graph $\mathcal{G}^{(\text{WM})}$ as a graph whose
vertex set $\mathcal{V}^{(\text{WM})}$ consists of all the white matter voxels
in a brain volume (approximately $250$k for the data used). The set of edges
$\mathcal{E}^{(\text{WM})}$ is defined such that each voxel is connected to
every voxel in a specified neighborhood, which we define as a $5\times 5\times
5$ region, resulting in at most $124$ neighbors for each white matter voxel.

We define the weight of the edge that joins a pair of voxels $(i,j)$ in
analogy to Iturria-Medina et al. [16] and Sotiropoulos et al. [17]. Let
$O_{i}(\omega)$ denote the diffusion orientation distribution function (ODF)
associated to voxel $v_{i}$, with its coordinate center being the center of
voxel $v_{i}$. Let $\mathcal{N}_{i}$ denote the set of vertices in
$\mathcal{V}^{(\text{WM})}$ that are adjacent to vertex $i$; i.e.,
$\mathcal{N}_{i}:\\{k\in\mathcal{V}^{(\text{WM})}|(i,k)\in\mathcal{E}^{(\text{WM})}\\}$.
For any two vertices $i,j\in\mathcal{V}^{(\text{WM})}$, let $\vec{r}_{i,j}$
denote the vector pointing from the center of vertex $i$ to the center of
vertex $j$, and define

$\displaystyle
p(i,\vec{r}_{i,j})=\int_{\Omega_{i,j}}O_{i}^{n}(\omega)d\omega,$ (2)

where $n\in\mathbb{Z}^{+}$ is a desired power factor to sharpen the ODFs and
$\Omega_{i,j}$ denotes a solid angle of $4\pi/98$ around $\vec{r}_{i,j}$
subtended at the center of voxel $v_{i}$. It is desirable to sharpen the ODFs
since they generally manifest only slight variations between directions of
strong and weak diffusion. Using (2), the weight between vertices $i$ and $j$,
denoted $w_{i,j}$ is defined as

$w_{i,j}=\frac{p(i,\vec{r}_{i,j})}{C_{i}}+\frac{p(j,\vec{r}_{j,i})}{C_{j}},$
(3)

where $C_{k}=2\max_{l\in\mathcal{N}_{k}}p(k,\vec{r}_{k,l})$. The normalization
factor ensures having $w_{i,j}\in[0,1],\quad\forall
i\in\mathcal{V}^{(\text{WM})},\forall j\in\mathcal{N}_{i}$.

Using diffusion MRI, $O_{i}(\omega)$ is estimated at a discrete set of
directions $\\{\vec{r}_{k}\\}_{k=1}^{N_{o}}$ from the center of an ODF,
denoted $\\{O_{i,k}\\}_{k=1}^{N_{o}}$, which can be used to approximate
$p(i,\vec{r}_{i,j})$ as a sum

$\displaystyle
p(i,\vec{r}_{i,j})\approx\frac{4\pi}{N_{o}}\sum_{k\in\mathcal{D}_{i,j}}O_{i,k}^{n},$
(4)

where $\mathcal{D}_{i,j}:\\{k\>|\>\vec{r}_{k}\in\Omega_{i,j}\\}$. The
resulting edge weights amount to a measure of coherence in the directions of
diffusion at neighboring voxels, with high weights associated with highly
coherent diffusion and vice versa.

###  2.3 Graph filter definition

Conventional filters defined within the Euclidean domain, such as Gaussian
filters, encompass a shift invariant impulse response. In contrast, vertex
realizations of graph filters are _shift-variant_ , leading to a unique
spatial realization of the filter at each graph vertex. This property,
combined with the weighting scheme in the proposed graph definition, results
in graph filters adapted to the microstructure of white matter.

Due to the lack of shift-invariance, it is convenient, and conventional, to
specify graph filters in the graph spectral domain. In this work, we design
and leverage smoothing filters associated with a heat kernel spectral profile,
defined by

$\quad\mathcal{K}(\lambda)=e^{-\tau\lambda},\quad\forall\lambda\in[0,\lambda_{\text{max}}],$
(5)

where $\tau$ is a free parameter determining the spatial extent of the filter.
In particular, $\mathcal{K}(\lambda)\in L^{2}(\mathcal{G})$, where
$L^{2}(\mathcal{G})$ denotes the Hilbert space of all square-integrable
${\mathcal{K}:[0,\lambda_{\max}]\to\mathbb{R}^{+}}$. Although such filters are
roughly analogous in shape to the Gaussian filters typically used for fMRI
data analysis, there is no direct equivalence between them.

Fig. 1: (a) Main orientation of ODFs manifesting white matter fiber directions
at a coronal slice in the vicinity of the corpus callosum. (b)-(d) Comparison
of the contour profiles of three spatial filters, with different extents of
adaptation to local anatomy, localized at a point within the corpus callosum
marked with a red dot in (a). (b) An isotropic Gaussian filter adapts neither
to local tissue morphology nor the underlying white matter structure. (c) A
filter restricted to white matter adapts to the local tissue morphology, but
not to the underlying white matter structure. (d) A filter defined on the
proposed graph adapts both to the tissue morphology and the underlying white
matter structure.

###  2.4 Graph filtering

Implementing graph filtering as in (1) requires calculation of all the
eigenvectors of the Laplacian matrix $\\{\bm{\chi}_{l}\\}_{l=1}^{N_{v}}$,
which is practically infeasible for larger graphs such as those proposed here.
Instead, we use a fast approximation algorithm [11]. Let $\mathcal{P}\in
L^{2}(\mathcal{G})$ be a polynomial approximation of kernel
$\mathcal{K}(\lambda)$. For a graph signal $\bm{f}$, its filtering with kernel
$\mathcal{K}$ can be found using $\mathcal{P}$ as

$\tilde{\bm{c}}_{\mathcal{K}}=\sum_{l=1}^{N_{v}}\mathcal{P}(\lambda_{l})\hat{\bm{f}}[l]\bm{\chi}_{l}=\mathcal{P}(\mathcal{L})\sum_{l=1}^{N_{v}}\hat{\bm{f}}[l]\bm{\chi}_{l}=\mathcal{P}(\mathcal{L})\bm{f},$
(6)

where $\tilde{\bm{c}}_{\mathcal{K}}\in\ell^{2}(\mathcal{G})$ with
$\tilde{\bm{c}}_{\mathcal{K}}[m]=\left(F_{\bm{k}}\bm{f}\right)[m]$; in the
last equality we exploit:
$\mathcal{L}\bm{\chi}_{l}=\lambda_{l}\bm{\chi}_{l}\Rightarrow\mathcal{P}(\mathcal{L})\bm{\chi}_{l}=\mathcal{P}(\lambda_{l})\bm{\chi}_{l}$.

This approximation has the benefit that it does not require the explicit
calculation of the eigenvectors. Instead, a polynomial of the Laplacian matrix
is applied to the signal, which can be efficiently implemented with matrix-
vector multiplication. Similarly to [11], we leverage a truncated Chebyshev
polynomial expansion approximation of $\mathcal{K}(\lambda)$, as it has the
benefit of approximating a minimax polynomial, minimizing an upper bound on
the approximation error.

##  3 Data

The MRI data used in this work was acquired from the Human Connectome Project
[18]. Structural images, parcellations, and diffusion data from the “$100$
Unrelated Subjects” subset of the $1200$ subject release were used. Data from
five of the subjects was excluded from the analysis due to incomplete white
matter coverage of the diffusion data. Additionally, fMRI data from a motor
task from one subject was used.

To evaluate the performance of the proposed graph-based filtering method we
created a set of semi-synthetic phantoms, consisting of synthetic activation
patterns generated on the basis of tractography results from real diffusion
data. Tractography was performed using generalized Q-sampling imaging [19] in
DSI Studio [20]. Synthetic activation patterns were created from individual
streamlines by diffusing an activation along the length of the streamline,
starting from a random point.

Activation patterns from $100$ streamlines were combined to produce a single
ground truth activation pattern for each subject. An example of such an
activation pattern can be seen in Figure 2(a). Time-series phantoms for all
subjects were then created by convolving the ground truth activation pattern
with a regressor corresponding to a block fMRI paradigm. The resulting time-
series were contaminated with additive white Gaussian noise of $\sigma=1$.

##  4 Results

###  4.1 Spatial adaptivity of graph filters

Figure 1 illustrates the differences in spatial adaptivity from various filter
definitions. Isotropic Gaussian filters adapt neither to the local tissue
morphology nor the underlying white matter microstructure. They transcend
tissue boundaries, integrating signal components associated to white matter as
well as the adjacent gray matter/CSF. Moreover, the filters integrate signal
components from multiple unrelated white matter tracts.

By restricting the analysis to the white matter tissue, the filters do not
transcend tissue boundaries and thus only combine signal components from white
matter. However, the filters can not differentiate signal components from
unrelated white matter tracts.

Filters defined on the proposed white matter graphs adapt both to the local
tissue morphology and the underlying microstructure. Their spatial profile
closely matches the anisotropic diffusion orientation manifested by the ODFs.

###  4.2 Semi-synthetic phantom results

The $95$ semi-synthetic time-series phantoms were filtered using both
isotropic Gaussian filtering and the proposed graph filtering approach. The
resulting filtered volumes were subjected to a standard GLM-based activation
mapping using the SPM toolbox. Figure 2 shows a comparison of the results
obtained by using the various filtering methods. For small filter sizes,
isotropic Gaussian filtering is capable of detecting the subtle shapes of the
activations, but with reduced sensitivity compared to graph filtering. More
activations are detected with larger filter sizes, but at the cost of
diminished specificity. In contrast, the shapes of the activations detected
through graph filtering are consistent across filter sizes.

The t-maps from all phantoms were thresholded at multiple levels and compared
with the available ground truth activation patterns in order to produce ROC
curves, which were then averaged across subjects. Figure 3 shows the average
ROC curves for a variety of filter sizes. The best performance was achieved
with Gaussian filters of $\text{FWHM}=2$mm and graph filters of $\tau=1.4$.
The proposed filtering approach results in increased sensitivity and
specificity across all tested filter sizes. Importantly, the performance of
Gaussian filters of $\text{FWHM}>2$mm is worse than without the application of
any filtering.

Fig. 2: (a) Example synthetic activation pattern. Dot-shaped activations
extend linearly in the plane orthogonal to the image. (b)-(e) t-maps obtained
from analysis conducted using: (b)-(c) isotropic Gaussian filtering,
$\text{FWHM}=2$mm and $6$mm respectively; (d)-(e) graph filtering, $\tau=1.3$
and $3.3$ respectively. All t-maps thresholded at $t=5$ and overlaid on the T1
image of the corresponding subject. Fig. 3: Average ROC curves from $95$ semi-
synthetic phantoms. The proposed graph-filtering approach shows increased
sensitivity and specificity across filter sizes.

###  4.3 Task fMRI results

As a proof of concept, results from a single HCP subject performing a motor
task are presented; see Figure 4. A substantially greater extent of
activations was detected in the corpus callosum from fMRI data filtered using
the proposed approach than from data filtered with isotropic Gaussian
smoothing. Similarly to the results obtained on phantoms, the detection
performance for isotropic Gaussian filtering deteriorates with increasing
filter sizes as a result of the increased ratio of extraneous signal mixed
with the signal of interest. In contrast, the proposed filters show more
homogeneous performance across filter sizes, as they consistently adapt to the
shape of the underlying signal.

Fig. 4: Single subject activations in the corpus callosum from a left hand
motor task using: (a)-(c) isotropic Gaussian filtering, $\text{FWHM}=2$mm,
$4$mm and $6$mm respectively. (d)-(e) graph filtering, $\tau=1.4$, $2.2$ and
$3.3$ respectively. All t-maps thresholded at 5% FDR.

##  5 Conclusion

We have proposed a novel graph-based approach for filtering white matter fMRI
data. The design enables constructing shift-variant spatial filters that adapt
to the underlying white matter structure, which enable revealing fine-grained,
anisotropic activity patterns. Results on semi-synthetic data showed the
potential of the proposed approach to enable improved specificity and
sensitivity in white matter activity detection, compared to the use of
isotropic Gaussian filtering. Our future work will focus on testing the
proposed approach on a larger cohort of subjects, to probe white matter
activations across different functional loading tasks as well as to
investigate resting state BOLD signal fluctuations in white matter.

## References

  * [1] Erin L Mazerolle, Steven D Beyea, Jodie R Gawryluk, Kimberly D Brewer, Chris V Bowen, and Ryan CN D’Arcy,  “Confirming white matter fMRI activation in the corpus callosum: co-localization with DTI tractography,”  Neuroimage, vol. 50, no. 2, pp. 616–621, 2010.
  * [2] Mara Fabri, Gabriele Polonara, Giulia Mascioli, Ugo Salvolini, and Tullio Manzoni,  “Topographical organization of human corpus callosum: an fMRI mapping study,”  Brain research, vol. 1370, pp. 99–111, 2011.
  * [3] Jodie R Gawryluk, Erin L Mazerolle, Steven D Beyea, and Ryan CN D’Arcy,  “Functional MRI activation in white matter during the Symbol Digit Modalities Test,”  Frontiers in human neuroscience, vol. 8, pp. 589, 2014.
  * [4] Matthew J Courtemanche, Carolyn J Sparrey, Xiaowei Song, Alex MacKay, and Ryan CN D’Arcy,  “Detecting white matter activity using conventional 3 Tesla fMRI: An evaluation of standard field strength and hemodynamic response function,”  Neuroimage, vol. 169, pp. 145–150, 2018.
  * [5] Yali Huang, Stephen K Bailey, Peiguang Wang, Laurie E Cutting, John C Gore, and Zhaohua Ding,  “Voxel-wise detection of functional networks in white matter,”  Neuroimage, vol. 183, pp. 544–552, 2018.
  * [6] Jodie R Gawryluk, Erin L Mazerolle, and Ryan CN D’Arcy,  “Does functional MRI detect activation in white matter? A review of emerging evidence, issues, and future directions,”  Frontiers in neuroscience, vol. 8, pp. 239, 2014.
  * [7] Zhaohua Ding, Yali Huang, Stephen K Bailey, Yurui Gao, Laurie E Cutting, Baxter P Rogers, Allen T Newton, and John C Gore,  “Detection of synchronous brain activity in white matter tracts at rest and under functional loading,”  Proceedings of the National Academy of Sciences, vol. 115, no. 3, pp. 595–600, 2018.
  * [8] Jodie R Gawryluk, Kimberly D Brewer, Steven D Beyea, and Ryan CN D’Arcy,  “Optimizing the detection of white matter fMRI using asymmetric spin echo spiral,”  Neuroimage, vol. 45, no. 1, pp. 83–88, 2009.
  * [9] Zhaohua Ding, Allen T Newton, Ran Xu, Adam W Anderson, Victoria L Morgan, and John C Gore,  “Spatio-temporal correlation tensors reveal functional structure in human brain,”  PloS one, vol. 8, no. 12, pp. e82107, 2013.
  * [10] Xi Wu, Zhipeng Yang, Stephen K Bailey, Jiliu Zhou, Laurie E Cutting, John C Gore, and Zhaohua Ding,  “Functional connectivity and activity of white matter in somatosensory pathways under tactile stimulations,”  Neuroimage, vol. 152, pp. 371–380, 2017.
  * [11] David K Hammond, Pierre Vandergheynst, and Rémi Gribonval,  “Wavelets on graphs via spectral graph theory,”  Applied and Computational Harmonic Analysis, vol. 30, no. 2, pp. 129–150, 2011.
  * [12] David I Shuman, Benjamin Ricaud, and Pierre Vandergheynst,  “Vertex-frequency analysis on graphs,”  Applied and Computational Harmonic Analysis, vol. 40, no. 2, pp. 260–291, 2016.
  * [13] Hamid Behjat, Nora Leonardi, Leif Sörnmo, and Dimitri Van De Ville,  “Anatomically-adapted graph wavelets for improved group-level fMRI activation mapping,”  NeuroImage, vol. 123, pp. 185–199, 2015.
  * [14] Hamid Behjat, Ulrike Richter, Dimitri Van De Ville, and Leif Sörnmo,  “Signal-adapted tight frames on graphs,”  IEEE Transactions on Signal Processing, vol. 64, no. 22, pp. 6017–6029, 2016.
  * [15] Antonio Ortega, Pascal Frossard, Jelena Kovačević, José MF Moura, and Pierre Vandergheynst,  “Graph signal processing: Overview, challenges, and applications,”  Proceedings of the IEEE, vol. 106, no. 5, pp. 808–828, 2018.
  * [16] Yasser Iturria-Medina, Erick Jorge Canales-Rodríguez, Lester Melie-García, Pedro A Valdés-Hernández, Eduardo Martínez-Montes, Yasser Alemán-Gómez, and JM Sánchez-Bornot,  “Characterizing brain anatomical connections using diffusion weighted MRI and graph theory,”  Neuroimage, vol. 36, no. 3, pp. 645–660, 2007.
  * [17] Stamatios N Sotiropoulos, Li Bai, Paul S Morgan, Cris S Constantinescu, and Christopher R Tench,  “Brain tractography using Q-ball imaging and graph theory: Improved connectivities through fibre crossings via a model-based approach,”  Neuroimage, vol. 49, no. 3, pp. 2444–2456, 2010.
  * [18] David C Van Essen, Stephen M Smith, Deanna M Barch, Timothy EJ Behrens, Essa Yacoub, Kamil Ugurbil, Wu-Minn HCP Consortium, et al.,  “The WU-Minn human connectome project: an overview,”  Neuroimage, vol. 80, pp. 62–79, 2013.
  * [19] Fang-Cheng Yeh, Van Jay Wedeen, and Wen-Yih Isaac Tseng,  “Generalized $q$-sampling imaging,”  IEEE transactions on medical imaging, vol. 29, no. 9, pp. 1626–1635, 2010.
  * [20] Fang-Cheng Yeh,  “DSI Studio,” http://dsi-studio.labsolver.org/,  [Online; accessed ].