^{1}

^{1}

^{2}

^{3}

^{†}

^{1}

^{2}

^{3}

^{*}

^{†}

^{1}

^{2}

^{3}

Edited by: Yalin Wang, Arizona State University, United States

Reviewed by: Niharika Gajawelli, Children's Hospital of Los Angeles, United States; Pew-Thian Yap, University of North Carolina at Chapel Hill, United States

This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience

†These authors share senior authorship

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Coordinate invariance of physical laws is central in physics, it grants us the freedom to express observations in coordinate systems that provide computational convenience. In the context of medical imaging there are numerous examples where departing from Cartesian to curvilinear coordinates leads to ease of visualization and simplicity, such as spherical coordinates in the brain's cortex, or universal ventricular coordinates in the heart. In this work we introduce tools that enhance the use of existing diffusion tractography approaches to utilize arbitrary coordinates. To test our method we perform simulations that gauge tractography performance by calculating the specificity and sensitivity of tracts generated from curvilinear coordinates in comparison with those generated from Cartesian coordinates, and we find that curvilinear coordinates generally show improved sensitivity and specificity compared to Cartesian. Also, as an application of our method, we show how harmonic coordinates can be used to enhance tractography for the hippocampus.

Traditionally, analysis of diffusion MRI (dMRI) images is performed in the Cartesian coordinates that the data is acquired in. Cartesian coordinates are well-suited for problems in which there is planar symmetry or in situations where there is no apparent symmetry or preferred direction. For example, in a MRI scanner the main magnetic field has an approximate planar symmetry, and usually scanner coordinates are chosen such that the

As we move our focus to specific areas within the body, the preferred directions depart from those offered by a simple Cartesian coordinate system, because most anatomical structures have complex and curved geometries for which curvilinear coordinates are more suitable. A typical example of this is the neocortex, which is a highly complicated structure due to its curved gyri and sulci. Since the microstructure of the cortex is arranged in laminae there is a preferred radial direction perpendicular to the surface of the cortex. The use of curved coordinates goes beyond just the choice of preferred directions. As shown by Bok (

For an application of the methods developed in this study we will focus on the hippocampus, a heavily studied structure of the brain. It plays a central role in learning processes, memory, and spatial navigation. It is also involved in major disease states (Hampel et al.,

The focus of our work will be dMRI tractography in curvilinear coordinates. Tractography is the only method for studying structural connectivity

The goal of this work is to compare tractography performed in Cartesian coordinates with tractography performed in curvilinear coordinates whose tangents point along the direction of the underlying microscopic fibres. We are particularly interested in configurations with high curvature as these are the ones that arise in the hippocampus and the neocortex. To make this comparison we perform simulations and construct the diffusion signal generated by a family of fibres that traverse around a sharp bend (tangential fibres) and are intersected by orthogonal fibres (radial fibres). To keep our comparison simple, and make interpretation of results more straight-forward, our simulations are restricted to two dimensions and we use conformal curvilinear coordinates. Next, by seeding only the tangential region we are able to gauge the performance of both coordinates by measuring the sensitivity and specificity of the resulting tracts. We find that in almost all scenarios curvilinear coordinates show improved performance over Cartesian coordinates.

As an application of our method we show how tractography can be performed on the hippocampus with the harmonic coordinates used in DeKraker et al. (

Let ^{i} = (^{i} is a scalar function of ^{j} = (^{j}) is straightforward to transform to new coordinates, we have simply, ^{i}): = ^{i}(^{j})). However, dMRI data is vector valued, i.e., for each b-vector,

Note that this transformation means that we would, in general, have different b-vectors for each voxel in an image. This requirement can be incorporated into existing dMRI analysis pipelines by composing the Jacobian matrix to the gradient non-linear distortion corrections. These distortions tend to warp the image in a spatially dependent way and are corrected using a Jacobian matrix. Many tools allow the option to include this Jacobian matrix (usually named

where,

These mathematical considerations that maintain generality allow one to perform tractography in arbitrary coordinates with the insertion of a few simple steps to existing pipelines. Namely, given the existence of new coordinates in terms of scanner-given Cartesian ones, we have to perform three steps, (1) resampling the dMRI data on a grid of the new coordinates, (2) Jacobian matrix transform of the b-vectors (either we use,

For simplicity and ease of interpretation we restrict our diffusion simulations to two dimensions. One straightforward way to generate conformal coordinates in two dimensions is by the use of the complex plane, ℂ. Let (^{−1}:(

where ϕ^{−1} ∈ ℂ, ^{−1}) and ^{−1}). The grid generated from this mapping is shown in

This figure summarizes aspects of our simulation strategy. In _{1/2}. In _{Σ}, is the seeding area, the region _{T} is the ground truth area occupied by the tangential tracts, and _{R} is used to detect tracts that falsely traverse into the radial fibres. Notice that _{R} is half of the region occupied by the radial tracts, this is done to reduce errors introduced by discretization. Also note that _{R} is asymmetric, this choice is to reject false positives caused by purely radial tracts emanating from _{Σ}, the focus is intended to be on the bent area. In

Next we simulate a diffusion signal from these coordinates. Along with two eigenvalues, λ_{u} and λ_{v}, we use the unit tangent vectors of these coordinate curves, _{ij}, we can generate a signal, _{1/2} = (_{min}, 0) + _{max}, 0))/2 be the midpoint of the bend at _{1/2} where _{1/2} = ℜ(ϕ(_{1/2}, 0)). To have the signal to transition from dominant tangential eigenvectors for _{1/2} to dominant radial eigenvectors for _{1/2}, we use a logistic windowing function, _{1/2})). We generate a tangential signal, _{T}, by choosing λ_{u} ≫ λ_{v} and a radial signal, _{R}, by choosing λ_{u} ≪ λ_{v}. Our final signal is then given by,

In our simulations we vary the following parameters, (1) resolution, (2) curvature of the bend with the ^{2}) and b-vectors of the first shell in the HCP dataset (Van Essen et al., _{0}, is 1,000, the dominant eigenvalue is 0.01 mm^{2}/s and the second eigenvalue is 0.0001 mm^{2}/s. The tractography step size is a quarter of the voxel length as per EuDx recommendation (Garyfallidis,

For all the tractography performed in this work we use the Constant Solid Angle ODF Q-ball model (Aganj et al., _{Σ}, as in shown _{Σ}. In addition, when we vary the resolution of our simulations, the mask for the seeds, _{Σ}, is kept at the highest resolution, this ensures that the seeding rate is the same across changes in resolution. For resampling all the diffusion data onto the new coordinate grid we use radial basis function interpolation.

To compare the tracts from the simulation obtained from using Cartesian vs. conformal coordinates we measure the sensitivity and specificity of the tracts for both approaches. The seeding region, _{Σ}, is labeled in _{0} to _{1/2}, _{T}, is also labeled. We also have the radial region, _{R}, which extends radially from _{3/4}, to _{max}. The reason to choose this sub-region over the whole region containing radial fibres is to avoid biases in the specificity caused by the discritization of the underlying grid, which is used to calculate areas. Also note in _{R} is slightly asymmetric under a reflection across the _{Σ} is omitted. This region is omitted to keep the focus of sensitivity and specificity on the bent region; since we use a windowing function and overlap the two families, there are tracts that are purely radial that emanate from the seed region and interfere with the sensitivity and specificity measurements in the region of high curvature. We use the following formulas to calculate the sensitivity and specificity,

We also calculate Youden's J-statistic, _{Σ} ∩ _{T}) is the number of voxels a tract passes through such that these voxels exist in both _{Σ} and _{T}, and Area(_{T}), on the other hand, is the area of the whole tangential region. Since we vary the resolution in our simulations, we choose a voxel size of 0.2 mm to calculate the area to reduce biases introduced by discretization. Occupancy (tract density) of tracts passing a voxel in a particular region is kept binary regardless of the number of tracts that pass through.

As an application of our methods we perform tractography with harmonic coordinates on the hippocampus. We choose two random subjects (four hippocampi) from the Human Connectome Project (HCP) Young Adult 3T study, WU-Minn S1200 release (Van Essen et al., ^{2} with approximately an equal number of acquisitions on each non-zero b-shell and an echo spacing of 0.78 ms. The resolution is 1.25 mm with isotropic voxels. Manual segmentation of all hippocampi and their subfields is performed, which is followed by solving for three harmonic coordinates, _{0} is the source and _{1} is the sink). The boundaries of these coordinates are motivated neuroanatomically, and consist of structures which border the hippocampus at its topological edge. We choose Neumann boundary conditions on the rest of the hippocampus to ensure that the tangent vectors of the coordinates are orthogonal to the surface normals at the boundary of the hippocampus. When solving for the coordinates _{0} and _{1} respectively), and Neumann boundary conditions on the rest of the hippocampus. More details about this procedure can be found in DeKraker et al. (

This figure summarizes the sensitivity and specificity measurements from the simulations. In

In

This figure shows the tracts generated from the simulations. For brevity we divide the resolution (

In _{0} is the source and _{1} is the sink, and Neumann boundary conditions are used for the rest of the hippocampus. An analogous procedure is used to compute the remaining coordinates. In

In this work we explored how changing from Cartesian coordinates to curvilinear coordinates affects the performance of tractography. If our data was continuous there would be no effect of changing coordinates. But since the problem at hand is discrete we see considerable changes brought about by a change of coordinates, when the constant coordinate curves follow underlying anatomical fibres. The effect we see is that the curvilinear coordinate system indirectly informs existing tractography methods of the underlying fibres, thereby acting as an anatomical prior. Other studies have also taken similar approaches to tractography, Dong et al. (

One of the key shortcomings of tractography done in Cartesian coordinates is in regions of high curvature (Schilling et al.,

As an application of our approach we performed tractography on four hippocampi with harmonic coordinates (DeKraker et al.,

One limitation of our work is that currently we are not able to mix coordinate systems, we can either perform tractography in Cartesian coordinates solely or curvilinear coordinates solely. In future work we will enhance our approach to mix coordinates systems while maintaining continuity of the tracts (and other mathematical objects), and thus are not limited to specific regions in the brain. In addition to simulations, we also demonstrated our approach in the

In this work we have presented a novel approach to tractography that utilizes prior knowledge about fiber pathways in the form of curvilinear coordinates to enhance the sensitivity and specificity of tracts, especially in regions of high curvature and low resolution. In addition, the procedure by which we have implemented our method is general, i.e., we can take existing tractography algorithms and easily integrate our approach into them. We have shown the enhancements achieved via robust simulations where, naturally, the ground truth is known, thus laying down a solid foundation for future work regarding curvilinear coordinates and dMRI.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below:

The studies involving human participants were reviewed and approved by Western University Health Sciences Research Ethics Board. The patients/participants provided their written informed consent to participate in this study.

UH: conception, methodology, code implementation, and writing of manuscript. CB and AK: conception, methodology, and manuscript editing. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Data were provided [in part] by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.