MesoFT: Unifying Diffusion Modelling and Fiber Tracking

1. Introduction

Diffusion MRI (dMRI) has become an essential tool for noninvasive mapping of brain tissue [5]. A unique advantage of dMRI originates from the Diffusion length, a typical displacement of water molecules, being of a few μm in a clinical scan, which is commensurate with cell dimensions. Hence, in addition to the millimeter-level anatomical MRI resolution, dMRI is sensitive to the tissue structure on the mesoscopic scale — an intermediate length scale between the molecular level where the NMR signal originates, and the macroscopic imaging voxel size. Mesoscopic brain tissue modeling attempts to quantify cellular-level tissue organization [2, 1, 4, 8] in each voxel. However, determining the μm-scale parameters, such as axonal dimensions, water fraction and myelin thickness from a dMRI signal acquired at a 100–1000 times lower resolution, is a difficult problem riddled by unstable model fitting and sensitivity to noise. Fiber tractography, on the other hand, focuses on characterizing the structural connectome and inferring the interregional relationships of the human brain. It delineates white matter tracts based on the empirical anisotropy of the dMRI signal, and does not attempt to quantify the mesoscopic structure. Tractography algorithms come in many varieties, divided into deterministic streamline-based, probabilistic and global [7, 9] approaches. Recently, the concept of tractometry [3] was developed, which projects mesostructural properties on top of already tracked fiber bundles. This idea combines both fields above, however it is just a rather a retrospective ‘combination’ of outcomes of conceptually different approaches.

In contrast, this work tries to merge both fields into one framework such that both problems can benefit from each other. There are a few attempts that try similar approaches [6, 10], but from an different perspective. The proposed algorithm performs a global tracking and local modeling simultanously. Instead of fitting the parameters independently in a voxel-by-voxel manner we treat the problem as one global optimization problem. On the one hand, we are able to link voxels by their structural relationships. On the other hand, the presence of fibers is purely driven by the fact that they correctly explain the observed signal and no further assumptions about topology are made. To make the optimization of this difficult, non-convex problem tractable we propose an efficient approximation of the data likelihood. A full brain reconstruction of the human brain takes about 12 hours on a standard Desktop PC.

2. The Fiber Model

The fiber model is built of small segments $X\in \mathcal{X}$. Each segment contributes to the predicted MR-signal $M(\mathcal{X})$ with a small signal contribution. Each segment carries its individual diffusion parameters that define this contribution. The segments can connect and polymerize to form long chains, called fibers. The set of edges connecting the segments is denoted by ε. The complete model $\mathcal{F}=(\mathcal{X},\mathcal{E},v)$ consists of the set of segments, their edges between them and the volume fractions υ. Our mesoscopic model M(r, q) is composed of axially-symmetric Gaussian Diffusion signals of the form

from different white matter compartments (here t is a fixed diffusion time and the b-value [5] b =|q|²t). The parameter n denotes the bundle direction, D_∥ and D_⊥ denote axial and radial diffusivities. The signal model is composed of the sum of two such tensor models, where for one of those the perpendicular diffusion is zero, and an additional constant reflecting non-diffusing water molecules. So, the signal from the ith segment is:

where υ_r is the volume fraction of completely restricted water, υ_a the intra-axonal fraction and υ_e = (1 − υ_a − υ_r) the extra-axonal fraction. While empirically this model captures the signal from a straight fiber bundle quite well, it currently has a limitation of setting the intra- and extra-axonal diffusivities to be the same. In general, these parameters seem not to be equal. The total expected signal is composed of a sum over all segments: $M(r,q)={\sum }_{X_i\in \mathcal{X}}{w}_{i}I\left(r,r_i\right)M_i(q)$, where I is an indicator function giving contributions if r and r_i is in same voxel. Each segment carries 5 variables $X_i=\left(r_i,n_i,D_{\Vert }^i,D_{\perp }^i,w_i\right)$, the position, direction, the axial diffusivity along the fiber, the perpendicular diffusivity and its overall weight. Note that the volume fractions υ_r(r), υ_a(r) and υ_e(r) = 1 − υ_a(r) − υ_r(r) are not properties of the segment but of the position. This is, on the one hand, conceptually quite natural and, on the other hand, a way to avoid the ambiguities of the model. To increase the number of segments (to get a higher number of fibers) the voxels may be divided into subvoxels which all share the same signal.

The cost functional, or energies as called in the following, consist of two parts: the data likelihood and the prior that regularize the problem and control the connections between the segments. For optimization we use, like in [9], a simulated annealing approach. The idea is to simulate the Gibbs distribution $P(\mathcal{F})=\frac{1}{Z}\text{exp}\left(-\left(E_{\text{data}}(M(\mathcal{X},v))+E_{\text{prior}}(\mathcal{F})\right)/T\right)$ while lowering the temperature T. For lower temperature it gets more and more likely that we sample from minimum of the energy. The simulation principle is based on a Reversible Jump Monte Carlo Markov Chain (RJMCMC).

3. The Algorithm

As already stated the optimization of the proposed energies is accomplished by an RJMCMC-type algorithm together with a cooling process. The idea behind the RJMCMC-algorithm is to repeatedly make random distortions to the current state $\mathcal{F}$. The distortion, called ${\mathcal{F}}^{\prime }$, usually depends on the previous state and follows some distribution $P_{\text{ prop }}\left(\mathcal{F}\mapsto {\mathcal{F}}^{\prime }\right)$, which can be arbitrarily chosen by the algorithm designer. The only condition is that the reverse transition has to be possible, i.e. $P_{\text{ prop }}\left({\mathcal{F}}^{\prime }\mapsto \mathcal{F}\right)>0$. The algorithm needs usually a certain number of initial iterations such that the sequence of generated states follows the desired distribution and is in equilibrium. Once equilibrium is reached (which can be checked by statistics of the energy differences), the system is slowly cooled down. In the following we present the different proposals used in our implementation.

4. Experiments

We consider a 2-shell scheme at b-values of 1000 and 2000 acquired with 60 directions per shell. The in vivo diffusion measurement was acquired on a Siemens 3T TIM Trio using an SE EPI sequence, with a TE of 107 ms. A healthy male volunteer (aged 36) was scanned at an isotropic resolution of 2.5mm. Additionally, a T₁ data set was acquired which was segmented into white matter (WM), gray matter (GM), and CSF using SPM. White matter was thresholded at a probability of 0.5 to determine the area of reconstruction.

First, to understand the importance of our approximation we did a brute force search on a synthetic data. By sweeping through the 3-parameter space of D_∥, D_⊥ and υ_i we found that our approximation speeds up the likelihood computation by a factor of 20 compared to an ordinary implementation. To validate the accuracy of the approximation we simulated a simple crossing/bending configuration (see Figure 1) consisting of three bundles. The central crossing has a crossing angle of 50°. The phantom is simulated on 24 × 24 × 9 grid with a isotropic voxel size of 2mm. Each of the three bundles has the same axonal volume fraction of 0.4, extra-axonal fraction of 0.6 and different diffusion parameters (D_∥, D_⊥). Bundle a) has (1, 0.5), bundle b) (1.5, 0.5) and bundle c) (2, 1). Rician noise was added with σ = 0.05 corresponding to a SNR of 20. Figure 1 shows histograms of our tracking results: fitted diffusion parameters, volume fractions, weight parameters and the tractogram. The reconstructed tracts are shown in three different coloring, one by directions, one by parallel diffusion D_∥ and one by perpendicular diffusion D_⊥. One can observe that all parameters are nearly unbiased. While the intra axonal volume fraction υ_a shows a small underestimation, the υ_r fractions and the diffusion coefficients show a small overestimation. For the in vivo dataset a voxel was subdivided into 3³ = 27 subvoxels to get a sufficient number of segments/fibers. For this setting the running time of the complete tracking procedure took about 10 hours on a Intel I7 (16GB) with four threads in parallel. The reconstruction contains 1.5 million particles forming about 50000 fibers longer than 10 segments.

In Figure 2 we show the results: Parametric brain maps of the diffusion parameters (Fig.2a), first and second order statistics (Fig.2b) of all parameters including tortuosity t = D_∥, D_⊥. The diagonal of the plot matrix shows ordinary histograms, the off-diagonal plots joint histograms of all parameter pairs. We also show histograms of the w parameter, the predicted signal at b = 0 and number of segments per voxels. Further, we selected several tracts (Fig. 2c) and d)) by two ROIs, namely, Cingulum (CG), Arcquate Fascicle (AF), Cortical Spinal Tract (CST), left Optic Radiation (OR), Fronto Occipital Fascicle (IFO) and callosal fibers to the precentral gyrus (CC). Finally, Figure 2e) shows fibers sliced coronally and colored by D_∥.

5. Discussion and Conclusion

We proposed a novel algorithm that unifies tractography and mesoscopic modeling to simultaneously reconstruct the human brain fiber bundle network and derives fiber specific diffusion parameters. The in vivo experiments show that the derived parameters go in-line with the current literature [4]. However, for the first time, we provide whole brain maps of the parameters including crossing regions. For single fiber voxel populations (like the Corpus Callosum) the putative axonal volume fraction (υ_a) is in a range of about 40 to 50 percent, while D_∥ ≈ 2 and D_⊥ ≈ 1 which is similar to [4], where these parameters where derived via kurtosis imaging. The inferred parameters from multi fiber voxels differ, one can observe reduced υ_a and D_∥ while an increase in D_⊥. The source of the restricted fraction υ_r is not yet clear. There is definitely a certain amount caused by the Rician noise. The generated tract bundles show similar appearance like [9], but they additionally carry individual diffusion parameters.