PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Inf Process Med Imaging. Author manuscript; available in PMC 2017 September 22.
Published in final edited form as:
PMCID: PMC5609858
NIHMSID: NIHMS870081

Topographic Regularity for Tract Filtering in Brain Connectivity

Abstract

The preservation of the spatial relationships among axonal pathways has long been studied and known to be critical for many functions of the brain. Being a fundamental property of the brain connections, there is an intuitive understanding of topographic regularity in neuroscience but yet to be systematically explored in connectome imaging research. In this work, we propose a general mathematical model for topographic regularity of fiber bundles that is consistent with its neuroanatomical understanding. Our model is based on a novel group spectral graph analysis (GSGA) framework motivated by spectral graph theory and tensor decomposition. GSGA provides a common set of eigenvectors for the graphs formed by topographic proximity measures whose preservation along individual tracts in return is modeled as topographic regularity. To demonstrate the application of this novel measure of topographic regularity, we apply it to filter fiber tracts from connectome imaging. Using large-scale data from the Human Connectome Project (HCP), we show that our novel algorithm can achieve better performance than existing methods on the filtering of both individual bundles and whole brain tractograms.

1 Introduction

Tractography is an essential tool in studying human brain connectomes using diffusion MRI (dMRI). It has been successfully applied to reconstruct important anatomical bundles and connectivity matrices for graph-based assessment of brain connectivity [1]. On the other hand, it is also well-known that tractography tends to generate a large number of false positives and false negatives [2]. While various tract filtering and clustering methods have been proposed, they are usually driven by the geometry of fiber tracts rather than their biological relevance [3, 4]. In addition, these methods rely heavily on the fine tuning of distance measures between tracts and heuristic selection of multiple parameters to achieve desired performance. Motivated by the wide presence of the topography-preserving property of anatomical connections in mammalian brains [5], we propose in this work a novel computational framework for tract filtering by measuring the topographic regularity of fiber tracts. Our method only involves the local analysis of topographic regularity and is generally applicable to both the reconstruction of individual fiber bundles and whole brain processing.

Topographic regularity is a common property in many brain networks [5]. The retinotopic organization of axons in the visual system, which involves almost half of the human brain, is an excellent demonstration of this property [6]. The retinotopy in the optic radiation is well-known from post-mortem studies [7]. Recent functional MRI studies showed that retinotopy is also followed in various pathways related to visual functions [5]. The somatotopic organization of the motor and sensory pathways [8], and the tonotopic organization of the auditory pathways [9] are also well-known examples where the fiber pathways follow a topograph reserving trajectory when they transmit sensory inputs to the cortical regions. Using data from tracer injection studies on macaque brains, the topographic organization of cortical connections to the striatum were also shown recently [10]. These regularities in brain pathways have great synergy to the recent proposition about the grid structure of fiber pathways [11].

While the topographic organization of brain pathways is a general principle, it has yet to be thoroughly investigated in tractography-based connectivity studies. A novel tractography-algorithm based on this principle was developed recently to model the topography-preserving connectomes in human brains [12]. In this work, we tackle the problem from a different perspective. We assume the fiber tracts have been generated from an arbitrary tractography algorithm and then filter the tracts based on their topographic regularity. Using tools from tensor decomposition, we quantify how topography is preserved as we traverse along each fiber tract. In contrast to the clustering approach that parcellate tracts into bundles, our method evaluates each tract separately by studying its topographic relation to a small set of neighboring tracts. In our experiments, we compare with spectral clustering methods, and demonstrate that our method can better improve the retinotopy of visual pathways and remove outlier tracts in whole brain tractography results.

2 Methods

2.1 Mathematical definitions and topographic regularity

In neuroscience, topographic mapping is generally defined as :point-to-point or region to region axonal connections that preserve spatial relationship between neurons [13, 5]. This definition implies that while nearby neurons are connecting to other nearby neurons, the spatial relationships among their projection trajectories do not vary. In this paper, we mathematically model this neuroanatomical property and apply it to filter tracts obtained using dMRI based tractography (Fig. 1). Below we describe our framework and definitions needed to mathematically study topographic regularity within this context.

Fig. 1
An example of a topographically regular bundle with the mathematical notation used in this work

Definition 1 (Fiber tract)

A fiber tract, or streamline, denoted as C={xi3|0iL} is a finite set of evenly spaced points along a curve.

Definition 2 (Tractogram)

A tractogram, T, is a set of fiber tracts, T={C0,C1,C2,}.

Definition 3 (n-tract neighborhood)

Given a fiber tract in a tractogram, C0T, the n-tract neighborhood of C0, is the set of n spatially closest fiber tracts to C0, including itself, and denoted by NC0={C0,C1,C2,,Cn1}.

In this work, we define the distance from tract Ck to tract Cl as the 1-sided Hausdorff distance, that is:

d(Ck,Cl)=maxxikCkinfxjlClxikxjl
(1)

where xjkxjl is the Euclidean distance between xik and xjl.xikCk, is the ith point on the kth fiber tract.

Definition 4 (n-point neighborhood)

Given a fiber tract, C0, its n-tract neighborhood NC0={C0,C1,C2,,Cn1}, and a point xi0C0, the n-point neighborhood of xi0 is defined as the set Nxi0={xi0,xi1,xi2,,xin1} such that xik,k[1,2,,n1] is the spatially closest points to xi0 on tract Ck.

Based on the n-point neighborhood, we can define the proximity measure between every pair of points within an n-point neighborhood as follows.

Definition 5 (Topographic proximity measure)

Given an n-point neighborhood of C0 around a point xi0C0, we propose a measure of proximity between any pairs of {xik,xil}Nxi0 as follows:

ρ(xik,xil)=ed(xik,xjk)/σ2
(2)

where σ is a model parameter and d is the normalized Euclidean distance computed as:

d(xik,xjk)=xik,xilmax{xikxil}Nxi0xikxil.
(3)

We call ρ the topographic proximity measure. In addition, we can define an n × n matrix Exi0 with elements Exi0(k,l)=ρ(xik,xil), and Exi0 is called the topographic proximity matrix for the n-point neighborhood of xi0

Note that the distance measure Eq. (3) and the topographic proximity measure Eq. (2) are rotation and translation-invariant. Besides, because of the normalization, they are scale-invariant as well. For example, for the fibers in Figs. 2(a) and (b) their topographic proximity matrices along the bundle are all identical despite the linear and nonlinear scaling. The twisting fiber bundle shown in Fig. 2(b) will also yield identical topographic proximity measures along the bundle due to the rotation invariance property. These properties agree with the topographic regularity concept in neuroanatomy. Finally we can define topographic regularity as follows:

Fig. 2
Topographically regular bundles with invariant affinity graphs.

Definition 6 (Topographic regularity)

A fiber tract C is topographically regular if Δ(Ex)=0xC or equivalently TGR=1mean(Δ(Ex))=1xC. Here 0Δ()1 is called topographic variation measure and TGR is called topographic regularity measure.

Our definition above interprets the neuroanatomical definition by quantifying the “spatial relationship between neuronal fibers” and the “preservation” of such relationship. We will clarify the mathematic definition of topographic variation and topographic regularity in the following section.

2.2 Spectral graph analysis for modeling topographic regularity

To define the variation of the topographic proximity matrices along any fiber, we should first quantify the distance between any two of such matrices. A topographic proximity matrix can be viewed as an affinity matrix of a graph formed by points within an n-point neighborhood of certain point on a fiber tract. In the literature of graph matching [14], it is often argued that the distance between the graph affinity matrices can be captured by the distance between the eigenvalues of their affinity matrices. In the following, we will exploit this idea.

Suppose we have two topographic proximity matrices, EX and EY. The distance between them can be written as:

ΔS(EX,EY)=EXEYF=SXSYF,
(4)

where SX and SY are eigenvalue matrices of EX and EY,F is the Frobenius norm. The last equality holds true if and only if the two matrices share the same set of eigenvectors. Eigenvalues for graphs are also known as the graph spectrum, and the associated analysis resides in the area of spectral graph theory [15]. However, graphs with identical spectra, namely isospectral/cospectral graphs, can be structurally distinct from each other. We refer the interested readers to [15] for details on this topic. Accordingly, we turn to look at the eigenvectors for graph modeling.

Eigenvectors have been proven critical in modeling graph structures in graph visualization [16] and graph partitioning [17]. We further demonstrate the char acteristics of eigenvectors in modeling topographical regularity via an example illustrated in Fig. 3. In this example, we first generate a non-rigidly distorted yet topographically regular fiber bundle, as shown in the top middle of Fig. 3. For the point cloud at the two ends of the bundle we can define the topographic proximity matrices, which form two graphs as shown at the top left and top right of Fig. 3. The sorted eigenvalues and their associated eigenvectors to each graph are shown underneath the respective graphs. We can observe that the eigenvalues are different and difficult to compare, while there is some obvious analogy in the top eigenvectors (highlighted by the dashed boxes). In addition, we also found two pairs of identical eigenvectors from their eigenvector sets. We have highlighted the identical eigenvector pairs in red and blue.

Fig. 3
Spectral graph analysis for nonrigidly deformed topographically regular bundle. The eigenvectors are arranged according to the magnitude of eigenvalues. The most dominant eigenvectors are shown at top left and the sequence goes from left to right and ...

In Fig. 4, we show another example of spectral analysis for a simple topographically irregular bundle. We only present the ends points due to limitted space. Because of the presence of a noisy fiber, we can observe a slight structural difference in the two graphs. After spectral analysis of the corresponding topographic proximity matrices, we observe certain minor differences in the eigenvalues of the two graphs. In addition, we observe that the top eigenvectors are identical while the last few eigenvectors are drastically different (magnitudes of their peaks and valleys do not match). These observations motivate us to model the topographical regularity using eigenvectors.

Fig. 4
Spectral graph analysis for a simple topographically irregular bundle. Left and right panels correspond to the two end points of the bundle.

Intuitively, we consider the following matrix distance based on eigenvectors:

ΔUY(EX,EY)=infoff(SX)=0EXUYSXUYTF
(5)

where off(A) is the off diagonal part of A, UY is the matrix formed by eigenvectors of Ey, and ΔUY(EX,EY) measures how well the eigenvectors of EY can be used to decompose EX. This distance measure is not symmetric, i.e., ΔUY(EX,EY)ΔUX(EX,EY), and, hence, it is biased. One way to symmetrize the measure is to use the ΔUY(EX,EY)+ΔUX(EX,EY). This method requires computing Eq. (5) (L − 1) times for any point along a tract consisting of L points. Therefore, it does not scale well for measuring topographic regularity of fiber tracts in massive tractograms. Alternatively, we adopt the following model:

ΔU(EX,EY)=infoff(SX)=0EXUSXUTF+infoff(SY)=0EXUSYUTF
(6)

where U is a matrix formed by a set of common eigenvectors for EX and EY. This method requires computing Eq. (5) only once for a point. It will be (L − 1) times faster than the aforementioned exhaustive approach.

This formulation is closely related to simultaneous diagonalization [18] which is a special case of Tucker tensor decomposition [19]. These methods try to find the optimal eigenvectors and eigenvalues at the same time, while our idea is slightly different in that we are more interested in the common eigenvectors and we consider the projected eigenvalues from the common eigenvectors as a byproduct, which greatly simplifies the computation.

Given a tract C={x1,,xL}, we can collect the topographic proximity matrices for all xiC as a 3rd-order tensor E={Exi|1iL}. To extract the common eigenvectors for the tensor E, we adopt the following model:

u=argmaxu=iEiuF2,s.t.uTu=1
(7)

where we used Ei as a simplified notation for Exi.

The intuition behind the above model is that we try to find the unit vector which is most similar to all the column vectors in the tensor. Since the column vectors are the topographic proximity measures, the optimal unit vector shall capture the dominant topographic proximity representation.

By differentiating the model and applying Lagrange multiplier, the optimal solution of the above is defined by:

u=λu,=[E1TE1E2TE2ELTEL].
(8)

Algorithm 1

Topographic regularity modeling of tractogram.

An external file that holds a picture, illustration, etc.
Object name is nihms870081f1.jpg

We shall call Eq. (7) group spectral graph model (GSGM). The above formulation also generalizes to multiple orthogonal vectors, and the optimal solution U* is defined by:

U=ΛU,UTU=I
(9)

where I is identity matrix, Λ is the diagonal matrix from Lagrange multipliers and it becomes the eigenvalue matrix when the equality holds true. We shall call the optimal solution U* the common eigenvectors of the tensor E. And we shall call Eq. (9) group spectral graph analysis (GSGA).

With the common eigenvectors, we can solve the back-projected eigenvalues:

Si=argminSiEiUSiUTF2,s.t.off(Si)=0
(10)

where off(Si) is the off diagonal part of Si. Note that EiUSiUTF2=UTEiUSiF2=diag(UTEiUSi)+off(UTEiU)F2, where diag(·) extracts the diagonal elements from the matrix that it operates on. The last equality is a basic property of Frobenius norm, i.e. AF2=ij|aij|2=i|aii|2+ij|aij|2. Resultantly,

Si=diag(UTEiU)
(11)

Based on Eq. (6) and the above derivations, the topographic variation measure for the fiber tract C can be formulated as:

Δ(Ei)=EiUSiUTFEiF,i=1,,L
(12)

where the matrix distance is normalized for scale invariance. The topographic regularity measure can then be computed as:

TGR=11LiLΔ(Ei).
(13)

The more regular the fiber bundle is, the higher TRG shall be. However, an issue with this formulation is that this form uses the norm of the residual and the residual may contain noise. We alternatively use the following conceptually equivalent yet cleaner formulation.

TGR=1LiUSiUTFEiF=1LiSiFEiF
(14)

The last equality holds because USiUTF2=Tr(Si2) by definition. Note that EiF=SiF and Si is the true eigenvalue matrix of Ei. The above measure is essentially the ratio of back-projected eigenvalues against the true eigenvalues.

Regarding eigenvector selection we have the following observation. In Fig. 5, we find that the dominant eigenvectors for each graph in Fig. 3 are also the common eigenvectors from GSGA, while the eigenvectors associated to small eigenvalues in the three cases differ from each other. This means only the common eigenvectors associated to the small eigenvalues can capture the irregularity occurred in the fiber tracts. Therefore, we only use the common eigenvector corresponding to the smallest eigenvalue in the GSGA for the entire fiber. Our algorithm is summarized in Algorithm 1. The first part of the algorithm is GSGA and the second half is regularity ranking and thresholding. The regularity measure we defined in Eq. (14) can be used to rank all the fibers and we shall remove the lowest η × 100% fibers according to the ranking.

Fig. 5
GSGA for Fig. 4. The top left and top right are the eigenvectors for the orange and green graphs in Fig. 4. Common eigenvectors from GSGA for both of the graphs are shown at the bottom.

3 Experimental results

We conducted both qualitative and quantitative evaluations to test our technique. In all the experiments we used data from the Q1-Q3 release of the Human Connectome Project (HCP) [20]. For all the experiments we computed fiber orientation distributions (FODs) using the algorithm in [21] and we used the iFODl algorithm of MRtrix3 [22] for tractography.

3.1 Qualitative evaluation

Visual pathway filtering

We reconstructed the optic radiation using the lateral geniculate nucleus (LGN) as seed and the primary visual cortex (V1) as the include regions. We used step size=0.1mm, angle=6°, cutoff=0.025,number=25K as tractography parameters. Spurious fiber tracts were removed according to fiber pairwise distance. Results were obtained using our method and spectral clustering based tract filtering [23] are shown in Fig. 6. We can observe drastic improvements of our method compared to the conventional spectral clustering. We find that the parameter combination of K = 100, σ = 0.01 and η = 0.2 gives overall the most satisfactory results which we will use in our quantitative evaluations on large-scale HCP datasets.

Fig. 6
Example results for filtering visual pathway of subject 100307. a) is the original optical radiation between V1 and LGN. b) is the result from spectral clustering. c)-f) are results of our method with different parameter settings. K is the neighborhood ...

In addition to outlier tract removal, our method also preserves topographically regular fiber tracts. Representative results are shown in Fig. 7. Our method does not only effectively remove outlier tracts but more importantly it also preserves topographically regular tracts.

Fig. 7
Example results for optical radiation bundles for subject 106319 viewed near V1. The image on the left is the original bundle, the middle is the result of the conventional spectral clustering based filtering, and the right one is our result. The dashed ...

Whole brain tractogram filtering

We also tested the applicability of our method on the whole brain tractograms. To obtain the tractograms, we used step size=0.1mm, angle=6°, cutoff=0.05, number=100K. In these tests, we compared our method with Spherical-deconvolution Informed Filtering of Tractograms (SIFT) [24]. SIFT is a tractogram post-processing approach that aims to find the best tractogram subset that matches the underlying FOD field. We obtained the SIFT output using the tcksift command of MRtrix3. The results are shown in Figure 8. From areas highlighted by the arrows, we can see that our method is more effective in removing spurious tracts running orthogonal to fibers immersing the cortex, thus significantly improving the topographical regularity in the filtered tracts than the SIFT filtered tractograms.

Fig. 8
Tract filtering for whole brain tractography (viewed near frontal lobe) of subject 100307 (top), 100408 (middle) and 101915 (bottom). Each row, from left to right, shows the original tractography result, the result of SIFT filtering, and the result of ...

3.2 Quantitative evaluation

The retinotopic organization of the visual cortex refers to a point-to-point mapping of the retinal space onto the visual cortex [6]. In [25], an automated technique was developed that assigns each vertex in the visual cortex two coordinates: angle and eccentricity. These coordinates have one-to-one correspondences on the retinal space. By leveraging this anatomical information, Aydogan and Shi recently proposed a validation measure of retinotopic organization for visual pathway fiber bundles based on the eccentricity component [12]. They observed that the eccentricity values form a “U” shape function of the cross-section of the optical radiation projecting onto V1. Thus, the goodness-of-fit measures for the quadratic eccentricity regression, such as mean-squared-error (MSE) and coefficient of determination (R2), are natural measures of topographic regularity for visual pathway fiber bundles. In this experiment, we also adopt these two measures for comparison. For quantitative evaluation we used 215 subjects from the Q1–Q3 release of HCP. The results are shown in Fig. 9. It is observed that our method consistently yields better fits that is measured with MSE. There is no significant difference in fit quality in terms of R2.

Fig. 9
Quantitative comparison of our method against spectral clustering based filtering for retinotopic organization. The MSE from our method is uniformly smaller and the R2 values are comparable.

4 Conclusion

In this work, we proposed a mathematical model of topographical regularity for tractograms. Our idea is to establish topographic representation for each point on the fiber and then use its variation to quantify the whole fiber’s topographic regularity. Using spectral graph theory and tensor decomposition, we proposed a mathematical topographic regularity definition that is well consistent with its neuroscientific understanding. We applied our model to filter visual pathway and whole brain tractograms obtained using data from HCP subjects. Both qualitative and quantitative evaluations show promising results compared to other filtering approaches. We believe our contribution with this work is significant. We not only translated the well-known neuroanatomical information of topographic regularity to mathematical formulation but also provided a novel tool for filtering tractograms using this definition; which makes our filtering technique fundamentally different compared to other approaches in the literature.

References

1. Bullmore ET, Bassett DS. Brain Graphs: Graphical Models of the Human Brain Connectome. Annual Review of Clinical Psychology. 2011;7(1):113–140. [PubMed]
2. Thomas C, Ye FQ, Irfanoglu MO, Modi P, Saleem KS, Leopold DA, Pierpaoli C. Anatomical accuracy of brain connections derived from diffusion MRI tractography is inherently limited. Proceedings of the National Academy of Sciences. 2014 Nov;111(46):16574–16579. [PubMed]
3. O’Donnell LJ, Westin CF. Automatic tractography segmentation using a highdimensional white matter atlas. IEEE Transactions on Medical Imaging. 2007;26(11):1562–1575. [PubMed]
4. Aydogan DB, Shi Y. Track Filtering via Iterative Correction of TDI Topology. Medical Image Computing and Computer-Assisted Intervention. 2015:20–27. [PMC free article] [PubMed]
5. Patel GH, Kaplan DM, Snyder LH. Topographic organization in the brain: searching for general principles. Trends in Cognitive Sciences. 2014;18(7):351–363. [PMC free article] [PubMed]
6. Engel SA, Glover GH, Wandell BA. Retinotopic organization in human visual cortex and the spatial precision of functional MRI. Cerebral Cortex. 1997 Mar;7(2):181–192. [PubMed]
7. Ebeling U, Reulen HJ. Neurosurgical topography of the optic radiation in the temporal lobe. Acta Neurochirurgica. 1988;92(1):29–36. [PubMed]
8. Ruben J, Schwiemann J, Deuchert M, Meyer R, Krause T, Curio G, Villringer K, Kurth R, Villringer A. Somatotopic Organization of Human Secondary Somatosensory Cortex. Cerebral Cortex. 2001 May;11(5):463–473. [PubMed]
9. Morosan P, Rademacher J, Schleicher A, Amunts K, Schormann T, Zilles K. Human Primary Auditory Cortex: Cytoarchitectonic Subdivisions and Mapping into a Spatial Reference System. Neuroimage. 2001 Apr;13(4):684–701. [PubMed]
10. Lehman JF, Greenberg BD, McIntyre CC, Rasmussen SA, Haber SN. Rules ventral prefrontal cortical axons use to reach their targets: implications for diffusion tensor imaging tractography and deep brain stimulation for psychiatric illness. Journal of Neuroscience. 2011;31(28):10392–10402. [PMC free article] [PubMed]
11. Wedeen VJ, Rosene DL, Wang R, Dai G, Mortazavi F, Hagmann P, Kaas JH, Tseng WYI. The Geometric Structure of the Brain Fiber Pathways. Science. 2012 Mar;335(6076):1628–1634. [PMC free article] [PubMed]
12. Aydogan DB, Shi Y. Probabilistic tractography for topographically organized connectomes. Proc MICCAI, Part I. 2016:201–209. [PMC free article] [PubMed]
13. Thivierge JP, Marcus GF. The topographic brain: from neural connectivity to cognition. Trends in Neurosciences. 2007 Jun;30(6):251–259. [PubMed]
14. Umeyama S. An eigendecomposition approach to weighted graph matching problems. IEEE transactions on pattern analysis and machine intelligence. 1988;10(5):695–703.
15. Brouwer AE, Haemers WH. Spectra of graphs. Springer Science & Business Media; 2011.
16. Kruskal JB. Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis. Psychometrika. 1964;29(1):1–27.
17. Shi J, Malik J. Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2000 Aug;22(8):888–905.
18. Bunse-Gerstner A, Byers R, Mehrmann V. Numerical methods for simultaneous diagonalization. SIAM Journal on Matrix Analysis and Applications. 1993;14(4):927–949.
19. Tucker LR. Some mathematical notes on three-mode factor analysis. Psychometrika. 1966;31(3):279–311. [PubMed]
20. Essen DV, Ugurbil K, et al. The Human Connectome Project: A data acquisition perspective. NeuroImage. 2012;62(4):2222–2231. [PMC free article] [PubMed]
21. Tran G, Shi Y. Fiber orientation and compartment parameter estimation from multi-shell diffusion imaging. IEEE Transactions on Medical Imaging. 2015 Nov;34(11):2320–2332. [PMC free article] [PubMed]
22. Tournier JD, Calamante F, Connelly A. MRtrix: Diffusion tractography in crossing fiber regions. International Journal of Imaging Systems and Technology. 2012 Mar;22(1):53–66.
23. Kammen A, Law M, Tjan BS, Toga AW, Shi Y. Automated retinofugal visual pathway reconstruction with multi-shell HARDI and FOD-based analysis. NeuroImage. 2016;125:767–779. [PMC free article] [PubMed]
24. Smith RE, Tournier JD, Calamante F, Connelly A. Sift: Spherical-deconvolution informed filtering of tractograms. NeuroImage. 2013;67:298–312. [PubMed]
25. Benson NC, Butt OH, Datta R, Radoeva PD, Brainard DH, Aguirre GK. The retinotopic organization of striate cortex is well predicted by surface topology. Current biology: CB. 2012 Nov;22(21):2081–2085. [PMC free article] [PubMed]