Introduction

Diffusion MRI (dMRI) data quality has improved dramatically in recent years, allowing more accurate estimates of orientation distribution functions (ODFs). However, accurate ODFs do not ensure accurate tractography, as different arrangements of axon bundles give rise to similar ODFs. In voxels with multiple ODF peaks, tractography must choose which one to follow arbitrarily (Fig. 1). This ambiguity is referred to as the "crossing vs. kissing fiber" problem but also encountered, e.g., where fiber bundles split or merge. We propose a solution that does not rely on prior anatomical information or on assumptions of microstructural continuity along bundles.

Methods

Local connection matrix (LCM): Fig. 2 illustrates the concept in 2D. Each square is a voxel in dMRI space, and the curves inside it are axon bundles. Let f_i, i=1,…,4 indicate 4 faces of a (2D) voxel. We define the LCM of the voxel as a 4x4 matrix C={c_ij}, where c_ij is non-zero only if there is a connection between voxel faces f_i and f_j. As C is symmetric, we only consider its upper triangular part. The figure shows three example bundle configurations with similar ODFs but different LCMs. Data acquisition: We used a human brain sample containing the area of where fibers split into the external capsule and internal capsules (EC/IC). The challenge for tractography is to follow these two bundles, but without following paths that go from the EC to IC or vice versa (Fig. 3). We acquired dMRI data at 9.4T (3D-EPI, 514-direction DSI¹, bmax=40K, 0.25mm isotropic resolution) and polarization-sensitive optical coherence tomography (PS-OCT)^2,3 scans on 63 coronal slices (0.003/0.075mm in-plane/through-plane resolution). The latter yield measurements of optic axis orientations. Ground truth: We performed tractography on PS-OCT orientation maps downsampled to 0.01mm. We counted the streamlines that connected each pair of faces in each 0.25mm dMRI voxel. For this proof of concept, we measured in-plane (2D) PS-OCT orientations. Thus we computed 4x4 LCMs, counting connections between the left, right, superior, and inferior faces of dMRI voxels. We normalized LCM elements to sum to 1, representing the portion of a voxel's fibers that connected each pair of faces. Network architecture: We constructed a convolutional neural network (CNN) whose inputs were the 3D diffusion propagators (zero-padded to 16³) from a target voxel and its 6 nearest neighbors. A convolutional encoder with 4 convolution layers of 60 features is followed by a ReLu non-linearity, batch normalization, and max-pooling. The deepest layer is flattened and dense connections are added to a layer with 512 elements, followed by a dropout layer (fraction=0.2) and a dense layer with a softmax activation function. The output comprises the 10 unique elements of the 2D LCM at the target voxel. We trained the model on data from half of the available slices, and used it to predict LCMs in the remaining slices. Tractography: We reconstructed ODFs from the DSI data. We performed streamline tractography, choosing which ODF peak to follow at each step 1. by the conventional rule of minimizing the bending angle with respect to the previous step, and 2. by sampling from the LCM.

Results

For each approach, Fig. 4 shows all streamlines that went through the EC (yellow sphere) or the IC (green sphere). No other inclusion regions were used. When we used the conventional choice of the peak that minimized the bending angle, results depended strongly on the maximum angle threshold. When this was set to a stringent 45^o (Fig. 4a), the EC bundle was not reconstructed fully, and a few erroneous streamlines turned from the EC into the IC. When we relaxed the threshold to 75^o (Fig. 4b), we recovered the EC bundle, but at the cost of more erroneous connections traveling from the EC to the IC and vice versa. Additional erroneous streamlines were found with both angle thresholds. When we used LCMs to determine which peak to follow, we did not find the aforementioned false positives. Fig. 4c shows the best-case scenario, using the ground-truth LCMs from the PS-OCT data. This illustrates that having access to the LCMs can indeed ensure that tractography takes the correct turn at each step. Fig. 4d shows the realistic scenario, usign LCMs predicted from the dMRI data by the CNN. This illustrates that we can recover LCMs successfully from DSI data. Importantly, when using the LCMs, we did not enforce a maximum angle threshold, i.e., we allowed tractography to take as sharp turns as it wanted, trusting the LCM to guide it to the correct choice.

Discussion

This early prototype used 4x4 LCMs to aid ODF-based tractography but, in the future, higher-dimensional (e.g., icosahedral) LCMs could replace ODFs altogether. We used full propagators from DSI as the CNN inputs, but will investigate if microstructural maps derived from less demanding (q-shell) sampling schemes are sufficient. This was a proof-of-concept experiment, where we trained and tested on different voxels from the same white-matter area. A larger and more diverse training set will be needed to scale this prototype up to the whole brain. We plan to build such a training set by obtaining realistic axon configurations from PS-OCT scans in different white-matter regions, and augmenting them by simulating additional configurations.

Acknowledgements

This work was supported by the National Institute of Biomedical Imaging and Bioengineering (R01-EB021265). It was carried out at the Athinoula A. Martinos Center for Biomedical Imaging at the Massachusetts General Hospital, using resources provided by the Center for Functional Neuroimaging Technologies, P41-EB015896, a P41 Biotechnology Resource Grant, and by the NIH Shared Instrumentation Grant Program (S10RR025563, S10RR023401, S10RR019307, and S10RR023043).