Keywords

1 Introduction

Tractography is a widely used technique for studying brain connectomes with diffusion MRI (dMRI) and has provided many exciting results in brain imaging research [1]. The lack of rigorous validation for in vivo human brain studies, however, has long been a critical challenge to push tractography toward a quantitative tool [2, 3]. On the other hand, the regular topographic organization of many fiber systems in human brains provide a surprisingly untapped anatomical knowledge for the improvement and validation of tractography techniques. Some of the well-known examples include the retinotopic organization of the visual pathway [4], the somatotopic organization of the somatosensory pathway [5], and the tonotopic organization of the auditory pathway [6]. In this paper, we incorporate this insight on anatomical regularity to develop a novel probabilistic tractography algorithm for studying the connectome of these topographically organized systems.

Conventional streamline methods rely on step size and curvature parameters to control the regularity of fiber tracks [7]. This type of local mechanism is difficult to account for the topographic organization of fiber pathways. It becomes particularly challenging when we want to reconstruct fiber tracks that are topographically arranged but with highly bended segments such as the Meyer’s loop of the optic radiation. More recently, global tractography techniques were developed to improve the robustness of streamline techniques [8, 9], but their focus is on balancing the fitting of dMRI signals and the regularity of local stick models. These developments are followed by microstructure based optimization techniques that start with a large set of tracks and assign weights or prune them so that resultant set of tracks best fit the dMRI signal [10, 11].

In this work we propose a novel probabilistic tractography method that incorporates both topographic and geometric regularity. Each fiber track is represented with the Frenet-Serret frame in our method. Fiber orientation distributions (FODs) [7, 12] are used in our work to represent the connectivity information at each voxel. For topographic regularity, the key idea in our method is the development of a novel likelihood function that tests the match of parallel curves to the FODs in the neighborhood of each point. The geometric prior is modeled as a Gaussian distribution of the Frenet-Serret frame. Taken together, we use a Bayesian approach with rejection sampling to propagate the fiber tracks and reconstruct highly organized fiber bundles.

In our experimental results, we validate the proposed technique on the reconstruction of the optic radiation using the multi-shell diffusion imaging data of the Human Connectome Project (HCP) [13]. The retinotopy of the visual system means there is a corresponding point on the cortex for each point of the retinal space [4]. This correspondence between the topography of the retina and axonal projections provides a great opportunity for localized mapping of retina disease to visual pathway integrity. Another distinct aspect of this fiber system is the Meyer’s loop that bends sharply toward the anterior aspect before moving toward the visual cortex [14]. Given these anatomical knowledge about the topography and geometric organization of the optic radiation, we believe it is an ideal testbed for tractography algorithms. On multi-shell imaging data from 56 HCP subjects, we compare our method with three tractography algorithms of the MRtrix package [7]. We demonstrate that our method is able to generate highly organized tracks while capturing the challenging Meyer’s loop. Using regression analysis, we quantitatively demonstrate that our method performs better in preserving the retinotopy of the optic radiation.

2 Methods

Frenet-Serret Apparatus: We use the Frenet-Serret frame and formulas to represent fiber tracks as differentiable curves. Let \(\gamma (s)\) represent a non-degenerate curve parameterized by its arc length s. The Frenet-Serret frame can then be defined with three orthonormal vectors: tangent (T), normal (N) and binormal (B) that are given as \(T=d\gamma /ds\), \(N=(dT/ds)/\mid dT/ds \mid \) and \(B=T \times N\). The Frenet-Serret formulas express the derivatives of T, N and B in terms of themselves with a system of first order ODEs: \(dT/ds=\kappa N\), \(dN/ds=-\kappa T+\tau B\) and \(dB/ds=-\tau N\). Here \(\kappa \) and \(\tau \) are the curvature and torsion of the curve, respectively. Given the initial conditions, this system can be uniquely solved.

Curve Spaces: Curve parameters can be put in a vector form to represent the space of differentiable curves that pass through \(p \in \mathbb {R}^3\) with \(\mathbf c _p=[\mathcal {F},\kappa ,\tau ]\in \mathcal {C}\). Here \(\mathcal {F}\) denotes the Frenet-Serret frame and \(\mathcal {C}\) is a short notation for the space of curves, \(S^2 \times S^2 \times \mathbb {R}^+ \times \mathbb {R}\), where \(S^2\) is the unit sphere. We will use the notation \(\gamma _\mathbf{c _p}(s)\in \mathbb {R}^3\) to denote where \(\mathbf c _p\) traces in space and denote \(\mathcal {F}_\mathbf{c _p}(s) = \{ T_\mathbf{c _p}(s), N_\mathbf{c _p}(s), B_\mathbf{c _p}(s)\}\), \(\kappa _\mathbf{c _p}(s)\) and \(\tau _\mathbf{c _p}(s)\) as the parameters of the curve at s.

Parallel Curves: Let \(\mathbf c _p\) and \(\mathbf c _{p'}\) be two curves. If \(p'\) lies on the normal plane of \(\mathbf c _p\) and \(\mathbf c _p=\mathbf c _{p'}\) then we call \(\mathbf c _p\) and \(\mathbf c _{p'}\) parallel curves and denote them as \(\mathbf c _p \parallel \mathbf c _{p'}\). These are also known as offset curves in computer graphics and computer aided design.

Fiber Model: We model a fiber track as a finite length, arc length parameterized curve that starts from a seed point \(p^{t=0}\in \mathbb {R}^3\) at time step zero (\(t=0\)). In a nutshell, our algorithm initializes by estimating \(\mathbf c _{p^{t=0}}^{t=0}\) for \(s=0\) and \(\gamma _\mathbf{c _{p^0}^0}(0)=p^0\). It solves the Frenet-Serret ODE system and moves to \(t=1\) by taking a step of \(\varDelta s\) which gives \(s=\varDelta s\), \(\gamma _\mathbf{c _{p^0}^0}(\varDelta s)=p^1\) and \(\mathbf c _{p^1}^{1}= [ \mathcal {F}_\mathbf{c _{p^0}^0}(\varDelta s), \kappa _\mathbf{c _{p^0}^0}(\varDelta s), \tau _\mathbf{c _{p^0}^0}(\varDelta s)]\). We then compute the priors, likelihood and posterior using Bayesian inference. Lastly we pick the next curve, \(\mathbf c _{p^1}^{2}\), randomly by rejection sampling and iterate until a stopping condition is met. Our track is a train of points traced by the curve \(\{p^0,p^1,p^2,\cdots \}\). Figure 1 explains our notation and the propagation technique.

Fig. 1.
figure 1

(a) At \(t=0\), without prior information, we select a random curve among the red candidate curves based on their likelihood. The thicker the curve, the higher posterior probability it has. (b) By solving the Frenet-Serret ODE, we propagate by \(\varDelta s\). (c) At \(t=1\), we calculate the prior probability for candidate curves for a smooth transition from the previous curve \(\mathbf c _{p^1}^{1}\) shown in green. (d) Propagate to \(p^2\).

Full size image

Bayesian Inference: Given a curve at time t, \(\mathbf c _{p^t}^t\), and data \(\mathcal {D}\), we estimate the posterior probability for the next curve \(\mathbf c _{p^t}^{t+1}\) using Bayesian inference as follows:

$$\begin{aligned} \underbrace{p(\mathbf c _{p^t}^{t+1} \mid \mathcal {D}, \mathbf c _{p^t}^t)}_{posterior} = \frac{p(\mathcal {D} \mid \mathbf c _{p^t}^{t+1}) p(\mathbf c _{p^t}^{t+1} \mid \mathbf c _{p^t}^{t}) }{p(\mathcal {D} \mid \mathbf c _{p^t}^{t})} \propto \underbrace{p(\mathcal {D} \mid \mathbf c _{p^t}^{t+1})}_{likelihood} \underbrace{p(\mathbf c _{p^t}^{t+1} \mid \mathbf c _{p^t}^{t})}_{prior} \end{aligned}$$
(1)

Prior Estimation: In order to preserve smoothness between the transition of curves that involves the changes in \(\kappa \), \(\tau \) and 3D rotation of \(\mathcal {F}\), we use Gaussian distributions to define the geometric prior. The variances for curvature and torsion are denoted with \(\sigma _\kappa ^2\) and \(\sigma _\tau ^2\). 3D rotation is realized by rotating \(\mathcal {F}\) around T, N and B axes in order, with variances \(\sigma _T\), \(\sigma _N\) and \(\sigma _B\) respectively. The overall prior probability is expressed in Eq. 2:

(2)

The functions used for computing the prior probability are given in Eq. 3.

$$\begin{aligned} \begin{aligned} p(T^{t+1} \mid T^t,\sigma _{N}^2,\sigma _{B}^2)&= (2\pi \sigma _N\sigma _B)^{-1}e^{-acos^2(\langle T', T^t \rangle )/2\sigma _{N}^2 - acos^2(\langle T^{t+1}, T' \rangle )/2\sigma _{B}^2} \\ p(N^{t+1} \mid N^t,\sigma _{T}^2,\sigma _{B}^2)&= (2\pi \sigma _T\sigma _B)^{-1}e^{-acos^2(\langle N', N^t \rangle )/2\sigma _{T}^2 - acos^2(\langle N^{t+1}, N' \rangle )/2\sigma _{B}^2} \\ p(B^{t+1} \mid B^t,\sigma _{T}^2,\sigma _{N}^2)&= (2\pi \sigma _T\sigma _N)^{-1}e^{-acos^2(\langle B', B^t \rangle )/2\sigma _{T}^2 - acos^2(\langle B^{t+1}, B' \rangle )/2\sigma _{N}^2} \\ p(\kappa ^{t+1} \mid \kappa ^t,\sigma _\kappa ^2)&= (\sqrt{2\pi }\sigma _\kappa )^{-1}e^{-( \varPsi (\kappa ^{t+1})-\varPsi (\kappa ^t) ) ^2/(2\varPsi (\sigma _\kappa ^2))} \\ p(\tau ^{t+1} \mid \tau ^t,\sigma _\tau ^2)&= (\sqrt{2\pi }\sigma _\tau )^{-1}e^{-( \tau ^{t+1}-\tau ^t ) ^2/(2\sigma _\tau ^2)} \end{aligned} \end{aligned}$$
(3)

In Eq. 3, \(\langle ~,~ \rangle \) is the dot product, \(T'\), \(N'\) and \(B'\) are the intermediary rotations, and \(\varPsi (\kappa )= asin(\kappa )\) is used to linearize change in curvature.

Likelihood Estimation: We estimate the data support for the next curve using the parallel curve definition and fiber orientation distribution (FOD). The field of FODs over an image volume can be expressed as a spherical function \(\mathcal {D}(p,T): \mathbb {R}^3 \times S^2 \rightarrow \mathbb {R}\), where \(p \in \mathbb {R}^3\) is the position of a point in the dMRI image and \(T \in S^2\). We define our likelihood expression as follows:

$$\begin{aligned} p(\mathcal {D} \mid \mathbf c _{p^t}^{t+1})= \frac{1}{4/3 \pi r^3} \int _{\forall \mathbf c _{p'}\parallel \mathbf c _{p^t}^{t+1} } \int _{\mid p^t-\gamma _\mathbf{c _{p'}}(s) \mid \le r} \mathcal {D}(\gamma _\mathbf{c _{p'}}(s),T_\mathbf{c _{p'}}(s))dsd\mathbf c _{p'} \end{aligned}$$
(4)

In Eq. 4, r is the radius of a sphere centered at p, which is the integration domain. Constant terms in front of Eq. 4 are used for normalization. The likelihood expression computes the data support for a candidate curve \(\mathbf c _{p}^{t+1}\) by integrating contributions of FODs belonging to parallel curves.

Implementation Details: The double integral in Eq. 4 is intractable. In order to estimate the likelihood numerically, we discretize the integral domain and add the contributions of FODs for a set of points as illustrated in Fig. 2. In this case, we used an actual human dMRI and picked a point in corpus callosum. For better visualization, we show here 7 random points as black dots within a radius of a few voxels, but in practice we use 27 points. The likelihood of a curve is computed by estimating the parallel curves that pass through these points. We then compute the tangents and interpolate the FODs at the points. The final likelihood is obtained by adding and averaging the data support contributed by these points. In Fig. 2(d) and (e), the estimated likelihoods of two different sets of parallel curves are visualized, where the thickness of fibers are in proportion to their likelihood.

Fig. 2.
figure 2

(a) To estimate the likelihood of a curve, we randomly pick a number of points within the integration radius, r. (b) Parallel curves passing through the random points. (c) We compute the tangents of parallel curves for each point and obtain the average FOD. (d)(e) Estimated likelihoods are shown in proportion to the thickness of curves.

Full size image

3 Test Subjects and Experimental Setup

In our experiments we used multi-shell images from Q1 release of HCP that includes 74 subjects, 56 of which completed both T1 and dMRI scans. For all 56 subjects, we computed FODs following the algorithm in [12] for multi-shell data. The FODs from this method are represented by spherical harmonics (SPHARM) and fully compatible with MRtrix that we compare with. We focused on the reconstruction of optic radiation that connects the lateral geniculate nucleus (LGN) and primary visual cortex (V1). One salient feature of this bundle is the retinotopic organization. Following the method in [16], we automatically generate V1 ROI and its retinotopic map that assigns each vertex in V1 cortex two coordinates: angle and eccentricity. ROI for LGN was generated using the method proposed in [14]. Inferior bundle of the optic radiation has an unconventional trajectory which first courses anteriorly before it runs posteriorly towards the visual cortex forming the elusive Meyer’s loop. Because of the quantitative coordinates provided by the retinotopic map and the challenges in reconstructing the Meyer’s loop, we picked the optical radiation in our tests. We conducted qualitative and quantitative comparisons of our technique with algorithms in MRtrix3 [7], which includes two probabilistic tractography algorithms: iFOD2 [15], iFOD1 and a deterministic, SD_STREAM, approach. Table 1 shows the parameters used for each technique. For our technique, we typically use a very small step size together with small \(\sigma _N^2\) and \(\sigma _B^2\) values since we walk along curves. In order to capture the Meyer’s loop, for MRtrix3 algorithms we used a higher angular threshold than the default parameters. Cut-off values were adjusted so that all the techniques can capture the Meyer’s loop. A lower cut-off threshold and a much higher angle threshold are used for SD_STREAM to achieve this.

Table 1. Tractography parameters used for each technique. vs is voxel size, \(^{\circ }\) is degree.
Full size table

4 Results

Qualitative Evaluation: Reconstruction results of the left optic radiation of an HCP subject by our method and MRtrix algorithms are shown in Fig. 3. We can clearly see that our results are more desirable as they are able to successfully capture the Meyer’s loop while exhibiting highly organized trajectories. As the tracks approach the V1 cortex, we can see the probabilistic tractography results from iFOD2 and iFOD1 start to become topographically less organized.

Fig. 3.
figure 3

Qualitative comparison of our method with MRTrix algorithms on the reconstruction of a left optic radiation of an HCP subject.

Full size image

Quantitative Evaluation: In Fig. 4 we illustrate our approach for quantitative evaluation. Using the eccentricity of the V1 cortex, we can divide it into three parts: fovea (red), superior-peripheral (green) and inferior-peripheral (blue) regions as shown in Fig. 4(a), which also allow us to split the fiber bundle into three sub-bundles. By cutting the bundle with a coronal plane (black dashed lines in Fig. 4(a)), we can visualize the topographic organization of the fiber tracks from different methods as shown in Fig. 4(b)–(e). Because of the topographic organization of the fovea and peripheral bundles, the eccentricity of the fiber tracks and their coordinates on the cross-section should follow a U-shape relation as plotted in Fig. 4(f), where the black dots are the raw data and the colored points are the fitted value with quadratic regression. In Fig. 4(g)–(j), the eccentricity values of each bundle on the cross-section are visualized, where the color bar for eccentricity is shown on the right most image. To quantitatively assess this relation, we applied quadratic regression to model the relation between eccentricity and the cross-sectional coordinates of fiber bundles. We report both the mean square error (MSE) and coefficient of determination (R\(^2\)) to measure how well the fiber tracks preserve the retinotopy of the fiber bundle. For each technique, the mean value of these two measures from 56 HCP subjects are listed in Table 2, where we can see that our method achieves the best performance in both measures.

Fig. 4.
figure 4

Quantitative comparison of the bundles show in Fig. 3. Top row shows the labels for three sub-bundles of the optic radiation. Bottom row shows the eccentricity values and the quality of quadratic fit using MSE and R\(^2\). The low MSE and high R\(^2\) values obtained by the proposed technique corroborate the qualitative observation.

Full size image
Table 2. Quantitative evaluation of the retinotopic organization of the optic radiation bundle. Best results are in bold.
Full size table

5 Discussions and Conclusion

Our current implementation is non-optimal and written in MATLAB and C++. Typically it takes several hours to generate tracks on the orders of thousands which we are working to improve. Our technique uses more parameters as compared to the conventional algorithms, but they are geometrically intuitive and their effects are not very difficult to understand. For quantitative comparison, we took a different approach than in [2] which suggests counting valid, invalid and no connections. This is because we focus on the reconstruction of individual bundles. We also did not use synthetic phantoms or simulated tracks because of the availability of retinotopic maps for in vivo validation.

In summary, we developed a novel probabilistic tractography technique that aims to capture the topographic organization of fiber bundles. A key idea in our method is the use of parallel curves to examine the local fitting of fiber tracks to the underlying field of FODs. Using the retinotopic mapping on V1 cortex, we have conducted quantitative evaluations and demonstrated that our method is able to generate more organized fiber tracks that follows known anatomy of the visual system. For future work, we will conduct more extensive validations on the visual pathway, its connectivity maps and other bundles that also follow topographic organizations such as the auditory and somatosensory pathways.