Diffusion MRI (dMRI) has become a very important tool for understanding the living brain tissue. Tractography tries to characterize the structural connectome to understand the details of the interregional relationships of the human brain at the macroscopic level. As the dMRI signal is a convolution of the microstructural fiber response function (FRF) with the mesostructural fiber orientation distribution (FOD),tractography algorithms usually also have to make assumptions about the FRF to infer the FOD. Most algorithms usually fix the FRF a-priori by analyzing single fiber voxels. Novel algorithms (2,3) also estimate the FRF during tractography. However, this also includes the microstructural modeling of the FRF, sometimes by simple tensor models, or by more complex multicompartement models. In this work, we introduce a novel way to infer FOD information with a minimal amount of a-priori assumptions about the nature of the FRF. The only assumption is that the dMRI signal in a single voxel is a convolution of the FOD with an axially symmetric FRF, no matter how the FRF looks. The corresponding energy for the FOD is non-convex and difficult to optimize, therefore we adopt the global tractography framework proposed in (3), which is based on a RJMCMC algorithm and has the flexibility to optimize the proposed objective.

In spherical harmonics the convolution of FRF and FOD takes the form of a product $$$S^b_{l,m}=M^b_l\,p_{l,m}$$$, where $$$M$$$ is the microstructural model (the FRF), and $$$p$$$ is the FOD and $$$b$$$ represents the dependence on a radial variable in q-space.For a given FOD $$$p$$$ we can easily optimize the quadratic objective $$$\sum_{b,l,m} |S^b_{l,m}-M^b_l p_{l,m}|^2$$$ for $$$M$$$ in unconstrained least-squares sense.If we now reinsert the optimal $$$M$$$ into the original quadratic objective we get the following energy $$ E(p) = -\sum_{l=0}^{l_\text{max}} \frac{\sum_b |\sum_m (p_{l,m})^* S^b_{l,m}|^2 }{ \sum_m |p_{l,m}|^2 }$$where all constants are neglected. If we now assume that $$$p$$$ is a superposition of weighted spherical delta peaks $$$p_{l,m} = \sum_{i=1}^N w_i Y_{l,m}(n_i)$$$, where $$$w_i$$$ is the weight and $$$n_i$$$ is the direction of the $$$i$$$th fiber segment, then the energy can easily be integrated into the global tractography framework proposed in (1). A fiber segment carries three parameter $$$(x_i,n_i,w_i)$$$, position, orientation and weight. The signal within a voxel is generated by all segments within the voxel. The position of the segment within the voxel does not influencethe signal generation. The number of segments is controlled by an additional $$$L_0$$$ penalty. The $$$L_0$$$ penalty is tuned such that the number of segments within a voxel reflects the crossing configuration, so usually not more than three segments per voxel. To increase the number of segments, the acquisitions voxels are subdivided into subvoxels.

We applied the proposed algorithm on a numerical bending/crossing phantom as illustrated in Figure 1. The data was sensitized at a b-value of $$$1000\ s/mm^2$$$ with different number of directions ($$$n=10,20,30$$$) and a SNR level of $$$20$$$.For $$$n=10$$$ the parameter $$$l_\text{max}$$$ was set to $$$2$$$, for $$$n=20,30$$$ to $$$4$$$. Figure 1 shows detailed results for $$$n=30$$$. Different ROI-based fiber selections are shown, the fiber/segment density of the tractogram, the terminal density and tract orientation density. Figure 2 shows a comparison for different $$$n$$$. In particular, for $$$n=10$$$ also the ROI-based fiber selections are shown. In-vivo we applied the algorithm to a usual $$$n=60$$$ scheme also at a b-value of $$$1000\ s/mm^2$$$ and compared to the original approach in (1). Figure 3 shows the comparison in terms of segment/terminal/tract-orientation density, Figure 4 shows the comparison for several selected tracts. The phantom results indicate that the proposed energy shows a unique ability to cope with a very small amount of gradient directions. Of course, the good results can be partly explained by the prior inside the global tracking algorithm prefering straight fibers. The reconstructed fiber densities and tract orientation densities (Figure 1e,1g) show very nicely the ability of the energy to deduce the right number of segments inside a voxel. This solves the problem which was inherent to the original algorithm (1), which produced lot's of 'broken' fibers around crossing regions.This is also shown in-vivo in terms of the terminal density. The hotspots of fiber terminals around the crossing regions are significantly reduced. The color-coded segment density shows the corresponding behavior of higher density in the crossing regions. Note that the proposed energy is not only suitable in the global tracking context, but can also used for simple single voxel analysis and is applicable to multishell data.