Brain Connectivity Measures via Direct Sub-Finslerian Front Propagation on the 5D Sphere Bundle of Positions and Directions

Abstract

We propose a novel connectivity measure between brain regions using diffusion-weighted MRI. This connectivity measure is based on optimal sub-Finslerian geodesic front propagation on the 5D base manifold of (3D) positions and (2D) directions, the so-called sphere bundle. The advantage over spatial front propagations is that it prevents leakage at omnipresent crossings. Our optimal fronts on the sphere bundle are geodesically equidistant w.r.t. an asymmetric Finsler metric, and can be computed with existing anisotropic fast-marching methods. Comparisons to ground truth connectivities provided by the ISBI-HARDI challenge reveal promising results, both quantitatively and qualitatively. We also apply the connectivity measures to real data from the Human Connectome Project.

Jorg Portegies, Stephan Meesters and Remco Duits: Joint main authors.

1 Introduction

We propose a connectivity measure between regions of the brain, based on diffusion-weighted MRI (dMRI) data. It is based on a quasi¹-distance defined through shortest paths in the joint space of positions and directions. Such paths are geodesics with respect to a highly anisotropic Finsler function (also known as ‘Finsler metric’). Although they can be recovered [9] from the quasi-distance that we compute numerically via anisotropic Fast Marching (FM) [26, 27], we do not require explicit tractography or streamline counting.

Tractography methods are commonly included in dMRI processing pipelines to model white matter anatomy (fiber tracts) in the brain. Tractograms can be used to define a measure for axonal fiber connectivity in terms of fiber counts between any two regions in the brain. Here we advocate that a connectivity measure should be derived directly from (wavefront propagations on) the 5D sphere bundle represented by dMRI data, rather than by counting streamlines in a tractogram.

For geodesic tractography, a global optimization problem is formulated based on a functional inferred from the dMRI data. Several functionals have been proposed, both for DTI and HARDI [11, 12, 19, 24, 25, 30, 37]. By virtue of their global nature, functional methods are more robust to noise than streamline methods. However, two issues must be dealt with, viz. proper handling of crossings and computational efficiency. Current front propagation methods typically take much longer than non-geodesic streamline tractography methods (cf. [6, 40]).

To this end we consider geodesic wavefront propagation on the sphere bundle, i.e. the base manifold $\mathbb {R}^{3}\times S^2$ of positions and directions, using the fiber orientation distribution function (FOD) obtained by constrained spherical deconvolution (CSD). CSD typically perfoms well in HARDI-challenges [4] and is included in neuroimaging software packages like DIPY and MRtrix [40]. Our connectivity measure relies on the quasi-distance mathematically scrutinized by Duits et al. [9]. Its four main characteristics are as follows:

• i. It includes a penalization of both length and curvature. Long and/or highly curved curves yield higher distances.

• ii. This penalization is based on the FOD. Curves poorly supported by the FOD yield higher distances.

• iii. It is defined on the 5D sphere bundle. This avoids leaking fronts at crossings, cf. Fig. 1 for an illustration on the $(2\!+\!1)$D sphere bundle.

• iv. It is induced by an asymmetric Finsler function defined on the sphere bundle. This allows for strictly forward geodesic front propagation, recall Fig. 1.

Fig. 1

Conventional geodesic wavefront propagation in the image domain (in red) typically leaks at crossings, whereas wavefront propagation in the domain $\mathbb {R}^{d}\times S^{d-1}$ (in green) does not suffer from this complication. For simplicity we depict the case $d=2$ in this figure. For dMRI the same idea applies to the case $d=3$. A minimum intensity projection over orientation gives optimal fronts in the image. The cost for moving through the orange parts is lower than elsewhere, and is computed from the FOD. The ‘leakage problem’ is absent both for propagating symmetric sub-Riemannian spheres (left) and for propagation of asymmetric Finsler spheres (right)

Theoretical Novelty. Our geodesic front propagation method for connectivity fits in the general contextual image processing paradigm used for denoising, fiber enhancement, and sharpening [5, 8, 28, 32, 34, 35, 39, 43]. It is similar to the pioneering work by Péchaud et al. [31]. There are, however, three major differences with our approach:

• i. In their definition of geodesic length [31], the anisotropy of the metric depends on the FOD, with isotropic behavior for low FOD amplitude. This was necessary to split numerical computations into a slow and fast part based on fast marching, respectively Dijkstra’s algorithm. However, in regions with low FOD values, we want to impose the same ‘stiffness’ of the curve for correct front propagation. Therefore, we propose a single, direct and accurate method that only computes optimal geodesic wavefronts.

• ii. It has now become feasible to do the entire computation using anisotropic FM (thanks to work by Jean-Marie Mirebeau [27], as an efficient alternative [36] to existing PDE models [2, 33]). The high efficiency of anisotropic FM comes along with only a limited and reasonable loss in accuracy [9, 33, 36], due to a special design of anisotropic acute stencils [26] that guarantee causality.

• iii. Our connectivity measure for source and sink regions $\mathcal {S}_{1},\mathcal {S}_{2} \subset \mathbb {R}^{3} \times S^2$ exploits asymmetric (sub)-Finslerian quasi-distances instead of (sub)-Riemannian distances, symmetrized afterwards by switching the role of $\mathcal {S}_{1}$ and $\mathcal {S}_{2}$ and taking the minimum. The driving asymmetric Finsler function, as opposed to a symmetric Riemannian metric tensor, recall Fig. 1, only allows for forward propagation. As a result it avoids highly unnatural cusps and it improves stability w.r.t. directions in source and sink regions, cf. [9, Def. 2, Thms. 2 and 3, Fig. 15] for proofs, visualizations of cusps, and further details.

The novelty over existing geometric path optimization approaches is as follows:

• i. Our geodesic wavefront propagation extends both Riemannian [11] as well as Finslerian wavefront propagation [24, 25, 37] to the larger base manifold $\mathbb {M}\doteq \mathbb {R}^{3}\times S^{2}$, i.e. the sphere bundle. This avoids overlapping fronts [37] (recall Fig. 1) and accounts for curvature.

• ii. We consider asymmetric Finsler functions and quasi-distances to enforce a spatially forward propagation, cf. Fig. 1, right column. Our framework also includes highly anisotropic symmetric Riemannian metrics for bi-directional propagation, see the middle column in Fig. 1. It admits infinite spatial anisotropy, constraining spatial motion at $(\mathbf {x},\mathbf {n}) \in \mathbb {R}^{3}\times S^{2}$ to the vector field $\mathbf {n}\cdot \nabla _{\mathbb {R}^{3}}= \sum _{i=1}^3 n^{i}\partial _{x^i}$. The latter relies on asymmetric sub-Finslerian or symmetric sub-Riemannian geometry. For convergence results, see [9, Thm. 2].

Experimental Contributions. After introducing our novel connectivity measure (accounting for crossings, stability and optimality) and its efficient implementation in Sect. 2, we demonstrate its use on two datasets:

• i. The ISBI HARDI Challenge 2013 dataset [4] provides a synthetic image simulating a complex configuration of fiber bundles with ground truth connections. We show that in the majority of cases, our connectivity agrees with true connectivity.

• ii. In human brain data from the Human Connectome Project (HCP) we estimate structural connectivity between the anterior nucleus of the thalamus and various cortical and sub-cortical regions. Results are generally in line with anatomical expectations.

2 Theory

We build on optimal path theory [9], using CSD [41] to construct an FOD image. We denote the sphere bundle by ${\mathbb M}\doteq \mathbb {R}^3 \times S^2$, and the FOD scalar field by $f : \mathbb {R}^3 \times S^2 \rightarrow \mathbb {R}^+$.

2.1 Measuring Curvature-Penalized Path Length

To obtain the length of a path $\gamma : [0,1] \rightarrow {\mathbb M}$ we construct a Finsler function ${\mathcal F}_0 : T{\mathbb M}\rightarrow \mathbb {R}^+$ on the tangent bundle $T{\mathbb M}$ over ${\mathbb M}$. For $\mathbf {p}= (\mathbf {x},{\mathbf n}) \in {\mathbb M}$ and $\dot{\mathbf {p}} = (\dot{\mathbf {x}},\dot{\mathbf {n}}) \in T_{\mathbf {p}}{\mathbb M}$ it is defined by:

$$\begin{aligned} {\mathcal F}_0(\mathbf {p},\dot{\mathbf {p}})^2 \doteq {\left\{ \begin{array}{ll} {\mathcal C}^{2}(\mathbf {p}) (\xi ^2 |\dot{\mathbf {x}} \cdot {{\mathbf n}}|^2 + \Vert \dot{{\mathbf n}}\Vert ^2) &{} \text {if } \dot{\mathbf {x}} = \lambda {\mathbf n}\,\text {with}\,\lambda \ge 0 \\ +\infty &{} \text {otherwise.} \end{array}\right. } \end{aligned}$$

(1)

The cost function ${\mathcal C}: {\mathbb M}\rightarrow \mathbb {R}^+$ will be explained in more detail later. The parameter $\xi >0$ balances the relative costs of spatial translation and in-place rotation. The length of a Lipschitzian curve $\gamma :[0,1] \rightarrow \mathbb {M}$ then equals

$$\begin{aligned} L(\gamma ) \doteq \int _0^1 {\mathcal F}_0(\gamma (t),\dot{\gamma }(t))\;\mathrm {d}t, \end{aligned}$$

(2)

and the (quasi-)distance between two points $\mathbf {p},\mathbf {q}\in {\mathbb M}$ is defined as

$$\begin{aligned} d(\mathbf {p},\mathbf {q}) \doteq \inf \{ L(\gamma ) \;|\; \gamma (0) = \mathbf {p}, \gamma (1) = \mathbf {q}, \gamma \in \text {Lip}([0,1],\mathbb {M}) \}. \end{aligned}$$

(3)

The distance is a specific solution (see [9]) to the sub-Finslerian eikonal PDE:

$$\begin{aligned} \sqrt{ \Vert \nabla _{S^{d\!-\!1}} U(\mathbf {q})\Vert ^2 + \xi ^{-2} ((\mathbf {n} \cdot \nabla _{\mathbb {R}^{d}}U(\mathbf {q}))_+)^2 } = \mathcal {C}(\mathbf {q}) \text { and }U(\mathbf {p})=0, \end{aligned}$$

(4)

where $U(\mathbf {q})\doteq d(\mathbf {p},\mathbf {q})$, $\mathbf {p}$ fixed, and $(a)_+\doteq \max \{a,0\}$. This sub-Finslerian PDE-system is a limit ($\epsilon \downarrow 0$) of the anisotropic Finslerian eikonal PDE system:

$$\begin{aligned} \left\{ \begin{array}{l} \sqrt{ \Vert \nabla _{S^{d\!-\!1}} U(\mathbf {q})\Vert ^2 + \xi ^{-2} ((\mathbf {n} \cdot \nabla _{\mathbb {R}^{d}}U(\mathbf {q}))_+)^2 + \epsilon ^{2} \Vert \mathbf {n} \wedge \nabla _{\mathbb {R}^{d}}U(\mathbf {q})\Vert ^2} = \mathcal {C}(\mathbf {q}) \\ U(\mathbf {q})=0 \text { if }\mathbf {q} \in A\doteq \{\mathbf {p}\} \, . \end{array} \right. \end{aligned}$$

(5)

We will admit source regions $A\supset \{\mathbf {p}\}$ and fix the anisotropy parameter $\epsilon =0.1$ as justified in previous work [9, 36]. The cost function ${\mathcal C}$ above is a transformation of the FOD field f that reflects the following conditions:

• Moving through an isotropic region should be expensive. To achieve this, we use on each position an $\mathbb {L}_1$ normalization on the sphere, i.e.

  $$\begin{aligned} f_1({\mathbf y},{\mathbf n}) \doteq \frac{f({\mathbf y},{\mathbf n})}{\int _{S^2} f({\mathbf y},{\mathbf n}) \mathrm {d}\sigma ({\mathbf n})} \text { and } f_2({\mathbf y},{\mathbf n}) \doteq \frac{f_1({\mathbf y},{\mathbf n})}{||f_1||_\infty } \in [0,1]. \end{aligned}$$

• To further increase the cost for motion through isotropic regions, we define an additional position-dependent (orientation-independent) cost ${\mathcal C}_{\text{ iso }}$:

  $$\begin{aligned} {\mathcal C}_{\text{ iso }}({\mathbf y}) = {\left\{ \begin{array}{ll} M &{} \max \limits _{{\mathbf n}\in S^2} f_2({\mathbf y}, {\mathbf n}) \le m \in [0,1]\\ 1 &{} \text {otherwise}, \end{array}\right. } \end{aligned}$$

  with $M=5$, $m=0.4$. We optionally include a white matter mask in $\mathcal {C}_{\text{ iso }}(\mathbf {y})$.

• Finally, we use the following cost function:

  $$\begin{aligned} {\mathcal C}({\mathbf y},{\mathbf n}) = {\mathcal C}_{\text{ iso }}({\mathbf y}) \frac{1+\sigma }{1+\sigma f_2({\mathbf y},{\mathbf n})^p} \text { with }p>1, \sigma >0. \end{aligned}$$

  (6)

  Increase of p sharpens the cost-function. Increase of $\sigma >0$ favors alignment with the ODF data-term over curve-stiffness.

2.2 Connectivity Measures

Connectivity between Points. Distance does not give information about connectivity directly, as it scales with the physical length of a path. We therefore define the connectivity $\kappa$ between two points $\mathbf {p},\mathbf {q}\in {\mathbb M}$ by a scale invariant ratio:

$$\begin{aligned} \kappa (\mathbf {p},\mathbf {q}) \doteq \frac{\int _0^1 {\mathcal F}_{{\mathcal C}\equiv 1} (\gamma ^*(t),\dot{\gamma }^*(t))\; \mathrm {d}t}{d_{{\mathcal F}_{{\mathcal C}}}(\mathbf {p},\mathbf {q})}, \end{aligned}$$

(7)

with $\gamma ^*$ the shortest path w.r.t. ${\mathcal F}_{{\mathcal C}}$, with ${\mathcal C}$ as above. This brings out the effect of the data along a path irrespective of its (Euclidean) length. Note that $\kappa = 1$ if ${\mathcal C}= 1$ along the entire geodesic $\gamma ^*$, i.e., if $\gamma ^*$ is perfectly supported by the data.

Connectivity between Regions. We define the connectivity $0 \le K(A,B) \le 1$ between regions $A,B \subset {\mathbb M}$ as follows: $K(A,B)=\min \{k(B,A),k(A,B)\}$, with

$$\begin{aligned} k(A,B) \doteq \frac{1}{|B|} \sum _{\mathbf {b}\in B} \frac{\int _0^1 {\mathcal F}_{{\mathcal C}\equiv 1} (\gamma ^*(t),\dot{\gamma }^*(t))\; \mathrm {d}t }{\min \limits _{\mathbf {a}\in A} d(\mathbf {a},\mathbf {b})} + \frac{1}{|A|} \sum _{\mathbf {a}\in A} \frac{\int _0^1 {\mathcal F}_{{\mathcal C}\equiv 1} (\gamma ^{**}(t),\dot{\gamma }^{**}(t))\; \mathrm {d}t }{\min \limits _{\mathbf {b}\in B} d(\mathbf {b},\mathbf {a})}, \end{aligned}$$

(8)

with $\gamma ^*$ and $\gamma ^{**}$ the minimizers of the denominators.

2.3 Fast Marching to Compute Distances

Propagation of geodesically equidistant wavefronts starting from a source point $\mathbf {p}\in {\mathbb M}$ is governed by the eikonal Eq. (4). With a numerical method called anisotropic FM [26, 27] one computes the geodesic distance to all other points $\mathbf {q}\in {\mathbb M}$ in a single pass.

Along with the geodesic distance (including the FOD-based cost) the curve length (excluding the FOD-based cost) can be computed without having to compute the actual optimizing curves, cf. [3]. Therefore, for each of the two terms in the connectivity measure (8) we need to run the FM solver only once.

The singular nature of the sub-Finsler metric results in numerical problems for FM solvers, which is why we need to approximate it with a highly anisotropic genuine Finsler metric, i.e. we solve (5) instead of (4). In theory we have convergence [9, Thm. 2]. In practice, this approximation is very close [9, 36].

3 Experiments and Results

Results for the ISBI HARDI Challenge 2013 Dataset. The ISBI HARDI Challenge 2013 dataset is a digital phantom that was designed for the ISBI 2013 Reconstruction Challenge. It represents a spherical volume that contains 27 bundles, and simulates dMRI with 64 uniformly distributed gradient directions with a b-value of $3000\,\text{ s/mm }^2$.

Since there are 27 bundles, each with a head (H) and a tail (T), we have 54 points on the boundary of the sphere (where the centerlines of these bundles cross the boundary). We assume these points to be known, and we run the FM algorithm 54 times, each time with a different point as seed (source) point and all other points as end points. This gives us 2809 geodesics that we indicate with $\gamma _{i_X\rightarrow j_Y}$, $i,j \in \{1,\dots ,54\}$, $X,Y \in \{H,T\}$, where the cases $i = j \wedge X = Y$ are excluded. The results of the FM and geodesic tracking are shown in Fig. 2. We see that the geodesics follow the anatomy of the dataset.

Fig. 2

Visualization of the ISBI dataset from different viewpoints. The two left-most images show the full dataset, the other four images show a selection of the bundles. The black lines inside these bundles indicate the geodesics obtained with FM

Fig. 3

Boxplot showing the sample median (white bar), 25% and 75% (gray box), max. and min. of connectivity values between the seed point (x-axis) and end points. Red dots indicate connectivity between seed and true end point. Since in most cases the red dot lies well above the gray box our connectivity measure is fairly distinctive

We compute the connectivity between each seed point and all end points, see Boxplot of Fig. 3 for the results. For example, choosing point ‘$1_H$’ as seed point (leftmost sample in the figure), we obtain 53 connectivity values with all other end points. The median, lower and upper quartile and maximum and minimum are computed from these values and plotted in the figure. The value of geodesic $\gamma _{1_H\rightarrow 1_T}$ (true connection) is indicated with the red dot. Ideally, this has the highest connectivity value. If we rank the geodesics by their connectivity values, then in 21 out of 54 cases the true geodesic receives the highest rank. The average rank by connectivity is 3.5 (1 being optimal). However, the method performs poorly on the two rather challenging cases $\gamma _{3_H\rightarrow 3_T}$ and $\gamma _{3_T\rightarrow 3_H}$ (if we exclude these, the average rank becomes 2.7).

Results for the HCP Dataset. The computation of the newly developed structural connectivity measure was evaluated for the connections of the anterior nucleus of the thalamus (ANT), and various cortical and subcortical atlas regions of healthy volunteers. The rationale for choosing the ANT as seed region was that this cortical structure (within the Papez circuit) may play a central role for epilepsy patients in the propagation of seizure activity, cf. [10, 14, 22]. Clinical trials have demonstrated that the DBS treatment is effective in reducing the median seizure frequency by 56% [10]. In a first preliminary step towards being able to assess why half of the patients are non-responders to the treatment, the proposed brain connectivity measures are applied to estimate the structural connectivity between the ANT and various cortical and subcortical atlas regions.

The HCP provides neuroimaging data from healthy subjects (N = 10) containing T1-weighted and T2-weighted anatomical brain scans as a well as a diffusion-weighted imaging (DWI) scan acquired by a multi-shell HARDI protocol (3 shells, 128 directions). Details on the scanning protocol can be found in work by Van Essen et al. [42]. The standard data preprocessing pipeline of the HCP project was utilized in which the usual corrections are applied [13]. The Morel atlas provides a collection of segmented thalamic nuclei in MNI space [17, 23, 29]. To segment the ANT for each individual subject, the anterior, medial and ventral subnuclei of the ANT, provided in the Morel atlas, were merged and a non-linear transformation to subject space was done by use of FSL [20]. Furthermore, the Harvard-Oxford atlas [7] was employed for segmenting various cortical and subcortical targets and was non-linearly transformed to subject space by the same procedure as the Morel atlas. The multi-shell CSD model [21, 41] available in the MRtrix software package [40] was used to compute a fiber orientation density (FOD) function of each DWI scan. The FOD function is subsequently aligned to the anatomical T1-weighted MRI via linear coregistration in FSL. We used the FOD function $\mathbf {p}\mapsto f(\mathbf {p})$ to define the cost function $\mathbf {p}\mapsto {\mathcal C}(\mathbf {p})$ according to (6), with $(p,\sigma )=(3,20)$, and the available gray-matter map to define the factor $\mathcal {C}_{\text{ iso }}(\mathbf {y})$. The bending-stiffness parameter was set to $\xi =0.1$. The FM algorithm is initialized with lifted seeds, i.e. points in ${\mathbb M}=\mathbb {R}^3 \times S^2$, using the sh2peaks function in MRtrix [40] to obtain the directional arguments.

As the full symmetric connectivity measure is computationally demanding, requiring a FM pass for each included atlas region, structural connectivity is computed only using the asymmetric connectivity metric (see r.h.s. of (8)):

$$\begin{aligned} k(A,B) \doteq \frac{1}{|B|} \sum _{\mathbf {b}\in B} \frac{\int _0^1 {\mathcal F}_{{\mathcal C}\equiv 1} (\gamma ^*(t),\dot{\gamma }^*(t))\; \mathrm {d}t }{\min \limits _{\mathbf {a}\in A} d(\mathbf {a},\mathbf {b})}, \end{aligned}$$

(9)

Here A denotes the ANT seed region and B the target cortical or subcortical brain region. The results of the connectivity analysis applied to healthy volunteers is presented in Fig. 5. The geodesic distance map is plotted (using our Visualization Tool available at https://goo.gl/D5Q7dE) relative to the ANT seed regions and target atlas regions. The average connectivity value is listed for several cortical and subcortical regions, relative to the ANT as seed region, computed for the left hemisphere (Fig. 5, bottom).

From anatomical and connectivity studies it is known that the ANT has extensive connections to the cingulate cortex and hippocampus [15, 16, 18, 38], to visual and frontal cortices [18] as well as a relation to the medial temporal lobe [1]. Indeed, it can be observed that the cingulate gyrus and the calcarine sulcus have a high connectivity (Fig. 5, top left). Furthermore, the frontal medial, temporal and insular cortices appear to have moderate connectivity values, where the inferior frontal gyrus and accumbens show lowest values.

Fig. 4

We visualize the geodesics $\gamma _{i_H\rightarrow j_T}$ and $\gamma _{i_T\rightarrow j_H}$, with color indicating connectivity-based rank, green (red) corresponds to high (low) rank. Akin to Fig. 3, many of the ‘true’ geodesics receive a higher connectivity value than ‘false’ geodesics

Fig. 5

Average connectivity value listed for a selection of cortical and sub-cortical regions from the Harvard-Oxford atlas, relative to the ANT as seed region, computed for the left hemisphere of 10 healthy volunteers. The anatomical locations of several regions is shown for a projection of the minimal geodesic distance map (top-left) and the T1-weighted anatomical image (top-right)

4 Conclusion

We proposed a new connectivity measure based on asymmetric sub-Finslerian (or symmetric Riemannian) wavefront propagation on the sphere bundle. It is capable of handling crossings in dMRI/HARDI data, uses only front propagation of optimal geodesics (without spatial cusps [9]), and can be computed with an accurate, direct anisotropic FM method (without the operator splitting in [31]).

We tested our connectivity measures on the synthetic data of the ISBI HARDI Challenge 2013. Comparisons to ground truth connectivities reveal promising results, both quantitatively (Fig. 3 in this article), and qualitatively (Fig. 4). We also tested our connectivity measures on HCP human brain data, where the structural connectivity results from the ANT region to surrounding ROI’s are in line with anatomical expectations (Fig. 5). These findings are promising in support of deep-brain stimulation for epilepsy treatment. However, in contrast to the phantom study, ground truth connectivity values are absent. Therefore, although the results on the HCP-data seem to indicate feasibility, this application should be further investigated.