Keywords

1 Introduction

Lots of mathematical models for the local diffusion process have been developed during the past decades from the diffusion tensor imaging (DTI) model to high angular resolution diffusion imaging (HARDI) models, hybrid diffusion imaging (HYDI) models and multi-compartmental microstructural (\(\mu \)-dMRI) models. Such models allow to compute ODF fields within a propagation domain and can be exploited to reconstruct virtual connections by following the optimal path corresponding to the peaks of ODF from a seeding position, yielding the so popular tractography concept that remains the unique method available to probe in vivo the human brain structural connectivity [1].

The simplest approach uses deterministic streamlining and consists in computing iteratively the backwards and forwards positions of a pathway by integrating the orientation information provided by the local diffusion model. This technique suffers from the low signal to noise ratio (SNR) inherent to the dMRI signal and corrupting orientational patterns, giving birth to artifactual connections. HARDI and HYDI models (see [2] for a detailed review) partly answered those issues and the emergence of probabilistic approaches helped in better managing corrupted or poor orientational information by allowing to temporarily track in sub-optimal directions [3]. Unfortunately, such methods still compute connections independently which was the main reason for the introduction of global tractography computing all the fibers at once. At that time, the computational cost for such algorithms was exceeding the capability of computers which is not the case anymore, giving a renew of interest for such approaches [4,5,6]. While improving robustness to acquisition errors or to ill-posed fiber configurations (crossing, kissing, fanning...), global tactography produces tractograms that still contain too many false positives due to the lack of anatomical constraints used to prevent their creation. Proposals have been made to increase the number of anatomical priors in order to overcome such issues [7, 8]. In this work, we propose a global spin-glass-based tractography algorithm relying on the anatomy to monitor the spins direction when entering gyri by allowing them to potentially connect to the cortex.

This paper is organized as follows: in the Methods section, we present the global tractography framework and our contribution, then we describe the optimization process. In the Results and discussion section, we demonstrate the efficiency of our global framework on a synthetic phantom and we highlight the inference of the structural connectivity from an ex vivo sample showing how well the fibers entering a gyrus may connect to its walls by following sharp turns. Finally, we discuss our results and further developments before concluding.

2 Materials and Methods

Part of global tractography approaches are spin glass-based methods consisting in optimizing a set of spins, each representing a small portion of fiber with respect to the local orientational pattern given by the dMRI data, and to other prior knowledge used to regularize and correct the corresponding connections. A spin glass is an oriented particle with a position \(x_{s}\), an orientation \(\mathbf n _{s}\) and a length \(2l_{s}\). We define a set S of N spins contained within the closed domain of white matter fascicles \(\varOmega \): \(S = \left\{ { s_{i}(x_{i},\mathbf n _{i},2l_{i}) | x_{i} \in \varOmega } \right\} _{0 \le i \le N}\). As described in [5], one can picture the spin glass as a small cylinder of null radius and with a length of \(2l_{s}\) (its forward extremity being at the position \(x_{s} + l\mathbf n _{s}\), and the backward one at \(x_{s} - l\mathbf n _{s}\)). What we hereafter call a spin glass connection is thus the association of two spins and each corresponding extremity: \(c_{1,2} = \big ( \left\{ s_{1}, \alpha _{1} \right\} , \left\{ s_{2}, \alpha _{2} \right\} \big )\) where \(s_{1} \in S, s_{2} \in S\), \(\alpha _{1} \in \left\{ +,- \right\} \) and \(\alpha _{2} \in \left\{ +,- \right\} \) (+ stands for the forward extremity and − for the backward). The corresponding set of spin glass connections is called \(\mathcal {C}\) and the global configuration M is therefore simply defined as the association of the spins and their connections: \(M=\left\{ S,C\right\} \).

2.1 Non-generative Anatomically Constrained Spin Glass Framework

In this work, we developed a novel global tractography method relying on a non-generative approach that creates a dense set of spins within a robust propagation domain [9], relies on any local diffusion model and is constrained by further anatomical priors. As in [5], we used Bayes’s theorem to model the a-posteriori probability to be optimized, yielding a global energy model to minimize \(P(M|D)=\frac{1}{Z}exp\big [- E_{global}(M) \big ]\) where Z is the partition function and \(E_{global}=E_{ext}(M,D)+E_{int}(M)\). The external energy \(E_{ext}(M,D)\), attached to the data D, and the internal energy \(E_{int}(M)\), of the global configuration M are detailed in the next paragraphs.

2.2 The Data Attachment Energy

With the design of our model, the external energy, or data attachment energy, should be minimum when the spin glass direction matches the direction of maximum diffusivity stemming from the acquisition, thus assumed to follow the direction of lowest restriction (e.g. the direction of fibers). Following [4], the global external energy can be written as:

$$\begin{aligned} E_{ext}(M,D) = \omega _{ext} \sum \limits _{S} -log\bigg [ p\left( x_{s},\mathbf n _{s} \right) \bigg ] \end{aligned}$$
(1)

where \(p(x_{s},\mathbf n _{s})\) is the probability of finding a spin glass at the location \(x_{s}\) and in the direction \(\mathbf n _{s}\). The design of our framework is generic, allowing to define this probability from any local diffusion model of ODF or PDF inferred from the adequate q-space sampling.

2.3 The Curvature Energy

For the internal energy, we extended the definition given by [5] adding a further dependency to the distance to the pial surface for the connection likelihood and the spin glass length. This allows to relax the curvature constraint in order to enable the creation of fibers depicting sharp turns when entering the cortex, as detailed in the next section.

$$\begin{aligned} E_{int}(M) = \omega _{int} \sum \limits _{c_{i,j} \in \mathcal {C}} U_{int}(c_{i,j}) \end{aligned}$$
(2)

where \(\omega _{int}\) is the internal energy weight and \(U_{con}(c_{i,j})\) is the interaction potential between two connected spins \(s_{i} = (x_{i},\mathbf n _{i}, 2l_{i})\) and \(s_{j} = (x_{j},\mathbf n _{j}, 2l_{j})\) defined in Eq. 3:

$$\begin{aligned} U_{int}(c_{i,j}) = \frac{||x_{i}+\alpha _{i}l_{i}{} \mathbf n _{i}-\overline{x}||^{2}}{l_{i}^{2}}+\frac{||x_{j}+\alpha _{j}l_{j}{} \mathbf n _{j}-\overline{x}||^{2}}{l_{j}^{2}} - \frac{(L_{i} + L_{j})}{2} \end{aligned}$$
(3)

with \(L_{i}\) and \(L_{j}\) referring to the local connection likelihood and \(\overline{x} = \frac{x_{i} + x_{j}}{2}\), being the middle point.

2.4 Embedded Anatomical Priors

At standard imaging resolutions, and because tractography aims at building fibers with low curvature, reconstructed fibers in a gyrus are generally limited to those connecting directly to its end, missing the ones connecting to its walls. Solving this issue is a key to get reliable connectivity matrices and address connectomics adequately. A first attempt was proposed in [10] relying on the use of a mean curvature flow algorithm performed on the pial surface to initialize the construction of streamlines using a deterministic approach. Following this idea, we propose an alternative strategy taking advantage of the global framework to constrain the spin glass length and direction during the initialization and the whole optimization process. It relies of the knowledge of the distance between the spin glass position and the pial surface inferred from the anatomical scan. Let’s note d this distance; then, the spin glass length can be written using Eq. 4, where \(l_{init}=0.8\frac{min(\varDelta _{x},\varDelta _{y},\varDelta _{z})}{N_{s}^{1/3}}\) with \((\varDelta _{x},\varDelta _{y},\varDelta _{z})\) being the mask resolution and \(N_{s}\) the number of spins initially created per voxel, meaning that it decreases as the spin glass gets closer to the cortex.

$$\begin{aligned} l_{s}=l_{init}\left( 1-exp(-\frac{d}{l_{init}})\right) \end{aligned}$$
(4)

Equation 5 defined the spin glass direction which was chosen as a trade-off between the orientation stemming from the local diffusion model and the normal direction \(\mathbf n _{pial}\) at the vertex of the pial surface corresponding to the projection of the spin glass position onto it, using a blending factor \(\psi (d)\).

$$\begin{aligned} \mathbf n _{s} = \psi (d) . \mathbf n _{pial} + (-1)^\beta \big [ 1 - \psi (d) \big ] . \mathbf n _{diff} \end{aligned}$$
(5)

\(\beta \) is equal to 1 or 2 in order to assign to the diffusion direction an outward orientation from the brain center to the cortex. \(\mathbf n _{diff}\) is computed using a Gibb’s sampler built from the probabilities corresponding to the angular profile of the diffusion model limited to the cone of axis \(\mathbf n _{s}\) and aperture angle empirically tuned (thus preserving the possibility to follow sub-optimal direction). \(\psi \) defined in Eq. 6 allows to cancel the impact of \(\mathbf n _{pial}\) on the spin glass direction when the orientation dispersion index ODI (as defined in [11]) is low, thus promoting the diffusion direction in this case, and enabling the spin glass to turn as the ODI increases towards \(\mathbf n _{pial}\).

$$\begin{aligned} \psi (d)= {\left\{ \begin{array}{ll} \frac{1}{2} OD \left( cos \left( \frac{\pi d}{d_{max}} \right) + 1 \right) , &{} if \ 0 < d \le d_{max} \\ 0, &{} otherwise \end{array}\right. } \end{aligned}$$
(6)

\(d_{max}\) controls the turn sharpness and corresponds to the sub-cortical ribbon where sharp turns occur. Thus, the connection likelihood L defined in Eq. 3 is modified accordingly, not to penalize the internal energy during sharp turn creation.

$$\begin{aligned} L=L_{init} + A log\left( \frac{d}{l_{init}}\right) exp\left( -\frac{2d}{l_{init}}\right) \end{aligned}$$
(7)

where \(L_{init}\) is similar to the one used in [5] and A set empirically to 25.

2.5 The Optimization Process

The optimization process was performed using a Metropolis-Hasting algorithm with a specific modification scenario designed to evolve. The available proposals are: (1) creation of a connection consisting in linking spins to form a fiber portion; (2) creation of a spin glass and connection prolongating a string of connected spins ending prematurely to progressively reach the boundary of the domain; contrary to the birth proposal defined in [5] linked to the generative nature of their algorithm, the creation of spins is not randomly defined here and aims at preventing the existence of fibers ending within the propagation domain; (3) optimal motion of spin glass clique, the clique corresponds to the chosen spin glass plus a limited set of its connected neighbors. Contrary to the optimal motion proposal defined in [5] that can lead to fibers oscillations if the energy path is too short, the proposal consists in modifying the position and orientation of the spin glass within the clique that minimizes the resulting clique tortuosity; (4) random motion of spin glass as defined in [5]; (5) death of spin glass. The modification proposal acceptance is driven by the following Green’s ratio \(GR=exp(-\frac{\varDelta E_{int}}{T_{int}} - \frac{\varDelta E_{ext}}{T_{ext}})\).

3 Results and Discussion

3.1 Numerical Fiber Crossing Phantom Experiment

A numerical phantom was designed to mimic a \(90^{\circ }\) crossing of two homogeneous population of fibers (where the branches can be assimilated to gyri) using a Gaussian mixture of tensors to generate a standard dataset of DW volumes corresponding to a spherical single-shell sampling of the q-space at b = 1500 s/mm\(^2\) and along 60 uniformly distributed diffusion directions over the sphere. A surface corresponding to the boundary was generated to mimic the pial surface. A region corresponding to the center of the crossing was drawn to mimic the ventricles and a distance map from any point of the propagation domain to this region was computed. An ODI map was designed to assign various OD indices to the different surfaces of the branches (from 0 to 1). The current global fiber tracking method was applied to the data using the analytical Q-ball model [12] (spherical harmonics order of 6, a Laplace-Beltrami regularization factor of 0.006) and the following parameters: 27 seeds per voxel; Gibb’s sampler temperature 0.015; initial external temperature 0.015; initial/final internal temperature 0.1/0.001; \(\omega _{ext}=100\), \(\omega _{int}=2000\); spin glass density 1.0; curvature threshold \(15.0^{\circ }\); maximum distance to pial 1.0 mm.

Fig. 1.
figure 1

Global spin glass tractograms obtained: (left) without any anatomical prior, connections are essentially from/to the top of the branches; (middle) with the pial surface and deep region distance map to constrain fiber directions close to the surface but reconstructing putative false positives fibers, connections are sparse across the propagation domain; (right) with the addition of the ODI map couterbalancing the systematic connection to the walls of the branches when spins get closer to the surface, the higher the ODI value is close to the surface, the more likely the connections to the walls of the branches.

Full size image

Figure 1 depicts the tractograms obtained using our algorithm without anatomical priors (left), taking into account the pial surface and the distance to the deep region (middle), and adding the ODI map (right). A connectivity matrix was computed from the center of the crossing to each surface and the results (shown mapped onto the mask) to evaluate the level of connectivity to the cortex. As expected, without anatomical priors, the connectivity matrix is sparse corresponding to the situation where fibers only connect the extremities of the branches. When the pial surface and deep region distance map are added, the connectivity matrix is the densest, due to the systematic creation of connections to the walls that might partially correspond to false positives induced by the proposed model. To counterbalance this effect, the ODI allows to control the creation of sharp turns with respect to the underlying anatomy. The density of the connectivity matrix is therefore decreased and the creation of sharp turns becomes dependent on the orientation dispersion level.

3.2 Real Human Brain Postmortem Sample Experiment

A human brain ex-vivo sample corresponding to a part of the visual cortex was scanned on a preclinical 11.7T MRI system equipped with a 780 mT/m gradient set and using a 60 mm volume coil. The imaging protocol includes a T2-weighted spin echo sequence (150 \(\upmu \)m isotropic) and a 2D PGSE sequence used to acquire a multiple shell diffusion-weighted dataset (300 \(\upmu \)m isotropic, b = 2000/4000/6000 s/mm\(^2\) with 60 / 125 / 125 diffusion directions respectively). The SHORE (polar harmonics order 4, Laplace-Beltrami regularization factor 0.006) ODF map [13] and NODDI (including extracellular, intra-axonal and stationnary compartments) ODI map [11] were computed from the DW dataset, and the pial surface and the ventricular region were extracted from the anatomical scan. The proposed global tractogaphy method was applied to the data using the same parameters as before except for the 1 seed per voxel in this case.

Fig. 2.
figure 2

High resolution T2-weighted postmortem scan (left) and the associated global spin glass tractogram obtained with anatomical prior knowledge

Full size image

Figure 2 depicts the resulting tractogram obtained with anatomical prior including the knowledge of the pial surface, the distance to ventricles, and the local ODI map (right) at the level of a gyrus, as depicted on the T2-weighted high resolution anatomical scan (left) using the red box. As expected, the addition of anatomical priors enables the inference of fibers connecting to the walls of the gyrus if assessed by a high level of orientation dispersion informed by the acquired diffusion data. Any connectivity matrix that would be computed using this novel approach would consequently be less contaminated by the overweight of connections to the top of gyri, thus allowing to use connectomics approaches more confidently to investigate the human brain structural connectome.

4 Conclusion

In this work, we have demonstrated that a spin glass global tractography framework is particularly suitable to easily embed any anatomical prior knowledge able to better constrain the creation of numerical connections. By doing so, one can more efficiently prevent the creation of false positives and deal with the local fibers configurations. The method takes advantage of the capability of a global approach to create all fibers at once to get an optimal tractogram, and also to adequately track fibers within gyri by creating connections not only at the end of them but also at their walls depending on the local orientation dispersion index. The efficiency of the method was demonstrated on both synthetic and postmortem human brain datasets. Future work will consist in adding further anatomical priors, for instance, to better monitor the seeding procedure, using microstructural features such as the neurite density, or to penalize the creation of connections between cortical and deep structures knowing their likeliness. Benchmarking our method onto a large database will be the next step to evaluate quantitatively its efficiency compared to alternative methods.