\authorinfo

Further author information: (Send correspondence to Zhangxing Bian: [email protected])

FastCod: Fast Brain Connectivity in Diffusion Imaging

Zhangxing Bian Department of Electrical and Computer Engineering, Johns Hopkins University, Baltimore, MD 21218, USA Muhan Shao Department of Electrical and Computer Engineering, Johns Hopkins University, Baltimore, MD 21218, USA Jiachen Zhuo Department of Diagnostic Radiology and Nuclear Medicine, University of Maryland School of Medicine, Baltimore, MD 21201, USA
Rao P. Gullapalli Department of Diagnostic Radiology and Nuclear Medicine, University of Maryland School of Medicine, Baltimore, MD 21201, USA Aaron Carass Department of Electrical and Computer Engineering, Johns Hopkins University, Baltimore, MD 21218, USA Jerry L. Prince Department of Electrical and Computer Engineering, Johns Hopkins University, Baltimore, MD 21218, USA

Abstract

Connectivity information derived from diffusion-weighted magnetic resonance images (DW-MRIs) plays an important role in studying human subcortical gray matter structures. However, due to the $O(N^{2})$ complexity of computing the connectivity of each voxel to every other voxel (or multiple ROIs), the current practice of extracting connectivity information is highly inefficient. This makes the processing of high-resolution images and population-level analyses very computationally demanding. To address this issue, we propose a more efficient way to extract connectivity information; briefly, we consider two regions/voxels to be connected if a white matter fiber streamline passes through them—no matter where the streamline originates. We consider the thalamus parcellation task for demonstration purposes; our experiments show that our approach brings a 30 to 120 times speedup over traditional approaches with comparable qualitative parcellation results. We also demonstrate high-resolution connectivity features can be super-resolved from low-resolution DW-MRI in our framework. Together, these two innovations enable higher resolution connectivity analysis from DW-MRI. Our source code is availible at jasonbian97.github.io/fastcod.

keywords:

Tractography, brain connectivity, thalamus parcellation

1 INTRODUCTION

Diffusion-weighted magnetic resonance imaging (DW-MRI) non-invasively images the diffusion of water in the brain, representing white matter (WM) tracts by depicting the anisotropy of the underlying microstructure. Based on DW-MRI, connectivity information can be obtained through probabilistic tractography. Connectivity is fundamental in studying the functionality and integrity of both cortical and subcortical gray matter structures (GM) [1]. The human thalamus, a subcortical GM structure that comprises numerous nuclei, functions as a central relay station for the brain, facilitating communication among sensory, motor, and associative brain regions. Since the nuclei are known to be connected to specific regions of the cerebral cortex, connectivity features can be used to parcellate the thalamus into nuclei [2, 1, 3, 4, 5, 6].

Despite the importance of connectivity features, the efficiency of computing connectivity is rarely studied. Traditionally, seeds are randomly sampled from the region of interest (ROI), from which millions of streamlines are generated by probabilistic tractography. The connectivity strength between ROIs is considered to be the number of streamlines that start from one ROI and end in another. Similar to the connectivity defined for two ROIs, we can define connectivity between voxels. For illustration purposes, we take the thalamic nuclei parcellation task as our test case; however, we note that our approach is applicable to many connectivity-based tasks [7, 8, 9, 10, 11].

Our goal is to find the connectivity of thalamic voxels and cortical regions. The complexity is $O(KNM)$ , where $N$ is the number of thalamic voxels which grows exponentially with the resolution of the image, $M$ is the number of targeted cortical regions which is proportional to the granularity of the chosen cortical atlas, and $K$ is the number of seeds per thalamic voxel for generating streamlines. It usually takes hours to extract connectivity for one subject. The inefficiency of current methods complicates obtaining high-resolution connectivity information (large $N$ ) or leveraging finer cortical atlases (large $M$ ). Also, $K$ is usually empirically chosen to be large enough to produce robust and reliable connectivity features [3], which inevitably increases the computational burden.

In this work, we propose a simple and accurate way to calculate connectivity. We define two ROIs or voxels to be connected if a streamline—no matter where it originates or ends—passes through them. We validated the proposed method on the Human Connectome Project (HCP) dataset[12] and show the proposed work achieves significant speedup without compromising the results. The speedup (of $120\times$ ) is more evident with high-resolution images. We also show we can super-resolve high-resolution features from a low-resolution DW-MRI.

2 METHOD

The connectivity feature for each voxel in the ROI (e.g., thalamus) $\mathcal{V}=\{v_{1},v_{2},\cdots,v_{N}\}$ is defined as the strength of its connection to $M$ targeted (e.g., cortical) ROIs. The strength of connection can be approximated by the number of streamlines connecting the voxel itself and the targeted ROI. The connectivity feature can be represented by a matrix $\mathbf{C}$ where the $i$ -th row of $\mathbf{C}$ is a $M$ -dimensional vector representing the connectivity of voxel $v_{i}$ to regions $\mathcal{R}=\{R_{1},R_{2},\cdots,R_{M}\}$ . Figure 1 shows the high-level idea that a connection is formed as long as a voxel is passed through by a streamline that connects to the target region.

Refer to caption — Figure 1: (a) A toy example. (b) The difference between the two approaches in calculating connectivities, given the total number of streamlines is the same, i.e., two streamlines in this case.

We compare the traditional (Alg. 1) and proposed (Alg. 2) algorithms shown below, where key differences are marked in red.

Input: $\mathbf{C}=\mathbf{0}$ ; $\mathcal{V}$ ; $K$ , the number of streamlines per voxel; target regions $\mathcal{R}$ Output : Connectivity matrix $\mathbf{C}$ 1 foreach $v_{i}\in\mathcal{V}$ do 2 $k\leftarrow 0$ 3 while $k<K$ do 4 Randomly sample seed $s_{j}$ inside voxel $v_{i}$ 5 if streamline $l_{j}\leftarrow\mathtt{Trac}(s_{j})$ is generated: then /* Streamline $l_{j}=\{p_{j}^{1},\cdots,p_{j}^{T_{j}}$ } */ 6 $r_{i}\leftarrow\mathtt{loc2region}(p_{j}^{T_{i}},\mathcal{R})$ 7 $\mathbf{C}(v_{i},r_{j})\leftarrow\mathbf{C}(v_{i},r_{j})+1$ 8 $k\leftarrow k+1$ 9 10 11 Algorithm 1 Traditional Approach Input: $\mathbf{C}=\mathbf{0}$ ; $\mathcal{V}$ ; $K^{*}$ , the number of total streamlines; target regions $\mathcal{R}$ ; $k\leftarrow 0$ Output : Connectivity matrix $\mathbf{C}$ 1 while $k<K^{*}$ do 2 Randomly sample seed $s_{j}$ in region $\mathcal{V}$ 3 if streamline $l_{j}\leftarrow\mathtt{Trac}(s_{j})$ is generated: then /* Streamline $l_{j}=\{p_{j}^{1},\cdots,p_{j}^{T_{j}}$ } */ 4 $r_{i}\leftarrow\mathtt{loc2region}(p_{j}^{T_{i}},\mathcal{R})$ 5 $\widetilde{\mathcal{V}}\leftarrow\mathtt{passthrough}(l_{j},\mathcal{V})$ 6 $\mathbf{C}(\tilde{v},r_{j})\leftarrow\mathbf{C}(\tilde{v},r_{j})+1$ for $\forall\tilde{v}\in\widetilde{\mathcal{V}}$ 7 $k\leftarrow k+1$ 8 9 Algorithm 2 Proposed Algorithm

In Alg. 2, instead of designating the number of expected streamlines for each voxel in $\mathcal{V}$ , we specify a total number of streamlines $K^{*}$ and randomly sample seeds in the source region $\mathcal{V}$ . In this way, every position in the source region has an equal chance to be seeded, and the resulting connectivity matrix directly reflects the underlying relative connectivity strength. Aiming at a fixed number of streamlines per voxel risks doing far more trials with seeding in inherently loosely-connected locations, leading to an overestimation of the connectivity.

The function $\mathtt{Trac}(s)$ represents the tractography algorithm which generates a streamline $l$ from seed location $s$ . A streamline $l$ consists of a group of ordered points { $p^{1},\cdots,p^{T}$ } in 3D, where $T$ is the number of points on this streamline. The function $\mathtt{loc2region}(p,\mathcal{R})$ returns the label of the target region $r\in\mathcal{R}$ where the point $p$ is located. In both algorithms, this function is used to determine which target region includes the endpoint of the streamline. In Alg. 2 line 5, the function $\mathtt{passthrough}(l,\mathcal{V})$ returns all the voxels $\widetilde{\mathcal{V}}$ that are passed through by the streamline $l$ . Then at line 6, the corresponding entries of the connectivity matrix record the connections. Lines 5 and 6 can be efficiently executed within Numpy [13] through array operations.

The proposed method has a complexity $O(K^{*}+NM/C)$ , where $C$ represents the average number of voxels in the source region passed through by a streamline. With an increase in resolution, $C$ will increase, but not as fast as $N$ does. Imagine an extreme case where every streamline passes through all thalamic voxels (i.e., $C=N$ ), then the complexity reduces to $O(K^{*}+M)$ . Note that the traditional algorithm (Alg. 1) has an inherent inability (due to inefficiency) to scale up to a higher resolution because the number of voxels $N$ in the source region exponentially increases with a higher resolution image, and it will quickly become computationally prohibitive. The proposed paradigm alleviates the computational cost by maximally exploiting each streamline’s contribution to the connectivity information $\mathbf{C}$ .

3 Experiments

Ten HCP [12] subjects are used for evaluation. The DW-MRIs and T1-weighted images are processed following the standard preprocessing pipeline [14] for HCP, resulting in thalamus segmentation and cortical parcellation. MRtrix3[15] is used to perform the probabilistic tractography. We set $K=200$ and $K^{*}=100,000$ .

Pie-glyphs. The output of our algorithm is an $N\times M$ connectivity matrix representing the connectivity of $N$ thalamic voxels to $M$ cortical regions. One common way to visualize connectivity results is to display a single color at each voxel where the color represents the region with the highest connectivity as shown in Figures 3 and 4. To gain a more complete appreciation of the connectivity, we present “Pie-glyphs”, as shown in Figure 2. Pie-glyphs show a pie chart at the center of each voxel location indicating the relative connectivity strength of that voxel to $M$ ROIs. Pie-glyphs provide less ambiguous information caused by partial volume effects, especially for large voxels, which is usually the case in diffusion-wegithed MRI. They also permit a more accurate appreciation of connectivity in voxels whose connectivity to two or more ROIs are approximately equal. For simplicity, however, in subsequent visualizations we use “single color” visual representation.

Comparision with traditional method. We tested both algorithms on ten subjects from the HCP dataset. The proposed algorithm is, on average, 30/48/85/124 $\times$ faster than the traditional method on 2/1.75/1.5/1.25 mm isotropic resolution with comparable qualitative parcellation results. Figure 3 shows one example. Interestingly, our approach generally produces more compact and less noisy parcellations. One possible reason for this is that the traditional approach requires a large $K$ to be robust enough against the randomness of the probabilistic tractography. However, using a large $K$ makes it even more computationally expensive. The processing time of our method does not increase exponentially with the resolution as does the traditional method.

Super-resolving connectivity features. One can upsample a streamline by inserting additional points between two neighboring points on the streamline. In theory, streamlines have infinite resolution under the assumption of smoothness. We perform tractography in the same resolution where DW-MRI is acquired (which is usually the lower resolution), but upsampling the streamlines by inserting evenly-spaced points into neighboring points, and then computing the connectivity in the structural image grid (which is usually the higher resolution). One can also use a smaller step size in the tractography algorithm, however, we find it is less efficient and does not make much difference. Since many neuroimage downstream tasks are performed in structural space, this super-resolving ability bridges the gap and allows us to directly compute high-resolution connectivity features. In practice, the rule of thumb for choosing the upsampling factor is that the spacing between two neighboring points on the upsampled streamline needs to be less than half of the voxel size of the resolution where connectivity is expected to be computed. For example, when the connectivity feature needs to be computed on 1 mm isotropic, the streamlines need to be upsampled to have $\leq 0.5$ mm spacing between points. Figure 4 shows two examples where the tractography is performed at a lower resolution (“baseline” results) and the connectivity features are extracted at a higher resolution (“SR” result).

4 Discussion

Calculating connectivity information is fundamental in studying the functionality and integrity of both cortical and subcortical GM structures. In this paper, we proposed an efficient algorithm to extract the connectivity features and tested it with the thalamus. We showed the proposed method is 30 to 120 times faster than the traditional method without compromising results. We also demonstrated a novel visual representation “pie-glyph” for connectivity information.

Acknowledgements

This work was supported by the National Institutes of Health / National Institute of Neurological Disorders and Stroke under grant R01-NS105503 (PI: R.P. Gullapalli). We want to thank Samuel W. Remedios (Johns Hopkins Univ.) for helpful discussion.

References

[1] Draganski, B. et al., “Evidence for segregated and integrative connectivity patterns in the human basal ganglia,” Journal of Neuroscience 28(28), 7143–7152 (2008).
[2] Behrens, T. E. et al., “Non-invasive mapping of connections between human thalamus and cortex using diffusion imaging,” Nature Neuroscience 6(7), 750–757 (2003).
[3] Stough, J. V. et al., “Automatic method for thalamus parcellation using multi-modal feature classification,” in [MICCAI ], 169–176, Springer (2014).
[4] Glaister, J. et al., “Thalamus parcellation using multi-modal feature classification and thalamic nuclei priors,” in [Proceedings of SPIE Medical Imaging (SPIE-MI 2016), San Diego, CA, February 27 – March 3, 2016 ], 9784, 97843J–97843J–6 (2016).
[5] Iglesias, J. E. et al., “Joint inference on structural and diffusion MRI for sequence-adaptive Bayesian segmentation of thalamic nuclei with probabilistic atlases,” in [IPMI 2019 ], 11492, 767–779 (2019).
[6] Yan, C. et al., “Segmenting thalamic nuclei from manifold projections of multi-contrast MRI,” in [Proceedings of SPIE Medical Imaging (SPIE-MI 2023), San Diego, CA, February 19 – 23, 2023 ], (2023).
[7] Hagmann, P. et al., “Mapping human whole-brain structural networks with diffusion MRI,” PloS one 2(7), e597 (2007).
[8] Gong, G. et al., “Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography,” Cerebral cortex 19(3), 524–536 (2009).
[9] Ye, C. et al., “A Bayesian Approach to Distinguishing Interdigitated Muscles in the Tongue from Limited Diffusion Weighted Imaging,” in [Bayesian & Graphical Models for Biomedical Imaging (BAMBI 2014) ], Lecture Notes in Computer Science 8677, 13–24 (2014).
[10] Bisecco, A. et al., “Connectivity-based parcellation of the thalamus in multiple sclerosis and its implications for cognitive impairment: A multicenter study,” Human brain mapping 36(7), 2809–2825 (2015).
[11] Shao, M. et al., “Direct reconstruction of crossing muscle fibers in the human tongue using a deep neural network,” in [Computational Diffusion MRI ], 69–80 (2021).
[12] Van Essen, D. C. et al., “The WU-Minn human connectome project: an overview,” NeuroImage 80, 62–79 (2013).
[13] Harris, C. R. et al., “Array programming with NumPy,” Nature 585, 357–362 (Sept. 2020).
[14] Glasser, M. F. et al., “The minimal preprocessing pipelines for the Human Connectome Project,” NeuroImage 80, 105–124 (2013).
[15] Tournier, J.-D. et al., “MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation,” NeuroImage 202, 116137 (2019).
[16] Desikan, R. S. et al., “An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest,” NeuroImage 31(3), 968–980 (2006).