# ORGANIZATION OF THE WHITE MATTER ANATOMY IN THE HUMAN BRAIN

EDITED BY : Laurent Petit and Silvio Sarubbo PUBLISHED IN : Frontiers in Neuroanatomy

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-314-2 DOI 10.3389/978-2-88963-314-2

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

# ORGANIZATION OF THE WHITE MATTER ANATOMY IN THE HUMAN BRAIN

Topic Editors:

Laurent Petit, Centre National de la Recherche Scientifique (CNRS), France Silvio Sarubbo, Ospedale Santa Chiara, Italy

Citation: Petit, L., Sarubbo, S., eds. (2019). Organization of the White Matter Anatomy in the Human Brain. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-314-2

# Table of Contents


Yue Bao, Yong Wang, Wei Wang and Yibao Wang

*29 Evidence for Plastic Processes in Migraine With Aura: A Diffusion Weighted MRI Study*

Nikoletta Szabó, Péter Faragó, András Király, Dániel Veréb, Gergő Csete, Eszter Tóth, Krisztián Kocsis, Bálint Kincses, Bernadett Tuka, Árpád Párdutz, Délia Szok, János Tajti, László Vécsei and Zsigmond T. Kincses

*36 Human Thalamic-Prefrontal Peduncle Connectivity Revealed by Diffusion Spectrum Imaging Fiber Tracking*

Chuanqi Sun, Yibao Wang, Run Cui, Chong Wu, Xinguo Li, Yue Bao and Yong Wang

*47 Awakening Neuropsychiatric Research Into the Stria Medullaris: Development of a Diffusion-Weighted Imaging Tractography Protocol of This Key Limbic Structure*

Darren W. Roddy, Elena Roman, Shane Rooney, Sinaoife Andrews, Chloe Farrell, Kelly Doolin, Kirk J. Levins, Leonardo Tozzi, Paul Tierney, Denis Barry, Thomas Frodl, Veronica O'Keane and Erik O'Hanlon


Daniel Schmitz, Sascha E. A. Muenzing, Martin Schober, Nicole Schubert, Martina Minnerop, Thomas Lippert, Katrin Amunts and Markus Axer

*91 The Superoanterior Fasciculus (SAF): A Novel White Matter Pathway in the Human Brain?*

Szabolcs David, Anneriet M. Heemskerk, Francesco Corrivetti, Michel Thiebaut de Schotten, Silvio Sarubbo, Francesco Corsini, Alessandro De Benedictis, Laurent Petit, Max A. Viergever, Derek K. Jones, Emmanuel Mandonnet, Hubertus Axer, John Evans, Tomáš Paus and Alexander Leemans

*109 A Missing Connection: A Review of the Macrostructural Anatomy and Tractography of the Acoustic Radiation*

Chiara Maffei, Silvio Sarubbo and Jorge Jovicich

# Editorial: Organization of the White Matter Anatomy in the Human Brain

Silvio Sarubbo<sup>1</sup> \* and Laurent Petit <sup>2</sup>

<sup>1</sup> Division of Neurosurgery, Structural and Functional Connectivity Lab Project, Azienda Provinciale per i Servizi Sanitari, Trento, Italy, <sup>2</sup> Groupe d'Imagerie Neurofonctionnelle, Institut des Maladies Neurodégératives (IMN)-UMR5293-CNRS, CEA, Université de Bordeaux, Bordeaux, France

Keywords: white matter, structural anatomy, functional anatomy, brain processing, neurosurgery, tractography, functional MRI, direct electrical stimulation

**Editorial on the Research Topic**

#### **Organization of the White Matter Anatomy in the Human Brain**

Between nineteenth and twentieth centuries, neurosciences experienced the first sharing of experiences and competences between the world of brain anatomy and clinics. The improvements in the knowledge of human white matter (WM) anatomy provided the natural background to the structural definition of a wide spectrum of clinical syndromes. This "disconnection" experience was the first field of strict integration between the WM anatomical and clinical skills, and constituted the hard core for the development of the modern neurosciences over the last century (Catani and ffytche, 2005).

While the second half of twentieth century has seen the neurophysiology taking a front role in the definition of the physiological and physio-pathological processing of brain circuitries, the last decade has definitively brought neuroimaging into the world of neuroscience. The functional magnetic resonance imaging (fMRI) and diffusion-weighted MRI (DWI) tractography have successively opened a new era for a better understanding of functional and structural anatomy of the human brain (Le Bihan and Johansen-Berg, 2012; Smith et al., 2013).

In particular, DWI-based tractography was the first tool allowing the exploration of human WM in vivo with an unprecedented level of details, and it shed a new light in the knowledge of the brain anatomy that became, finally, more accessible (Jeurissen et al., 2019). Beyond the technical aspects related to the continuous necessary improvement of this approach (Maier-Hein et al., 2017), tractography produced a conceptual revolution leading that the wiring diagram of brain connections regained a center scene of neuroscience research. Such a revolution was not only in research but also in the clinical and neurosurgical domains and opened the "connectome" era (Sporns, 2013).

The fields of neuroanatomy, neuroimaging, neurophysiology and clinical researches are currently closer as never before. In fact, two decades of exploration of brain structure and functional processing with an unprecedented level of sensitivity opened new challenges. Among others, the research for a ground truth in structural anatomy is definitely the most impressive, especially considering the basic and conceptual consequences of that in assessing a reliable knowledge of brain processing, clinics and plasticity. This is what the vast majority of the articles in this Research Topic highlight by describing association WM pathways (Bao et al.; David et al.; Panesar et al.), cortico-striatal Cacciola et al. and cortico-thalamic (Maffei et al.; Roddy et al.; Sun et al.) projection pathways.

Many techniques are now available for the exploration of structural anatomy of human WM. At macrostructural level, besides the millimetric and widespread power of resolution of DWI, polarized light imaging (PLI, see in this Research Topic the report of Schmitz et al.) and

Edited and reviewed by:

Javier DeFelipe, Cajal Institute (CSIC), Spain

#### \*Correspondence:

Silvio Sarubbo silviosarubbo@gmail.com

Received: 27 August 2019 Accepted: 04 September 2019 Published: 18 September 2019

#### Citation:

Sarubbo S and Petit L (2019) Editorial: Organization of the White Matter Anatomy in the Human Brain. Front. Neuroanat. 13:85. doi: 10.3389/fnana.2019.00085

**4**

polarization sensitive optical coherence tomography (PSOCT) are emerging techniques with an unprecedented level of accuracy in defining at micrometer resolution the structural anatomy of WM in humans (Caspers and Axer, 2017). Even, microdissection of fixated specimens (i.e., modified Klingler's preparation, Sarubbo et al., 2015b) experienced a renewed interest in confirming previous evidences or re-defining course and terminations of associative, commissural and projection pathways in human and non-human primate brain (Sarubbo et al., 2013, 2016, 2019; Fernandez-Miranda et al., 2015; De Benedictis et al., 2016; Hau et al., 2016, 2017; Wang et al., 2016; Maffei et al., 2018). Limitations and strength points of each of these techniques allowed to partially reduce the technical gap for clarifying and distinguishing among the open questions in the field of human WM neuroanatomy. Nevertheless, no single technique can be considered alone as fully reliable for a whole and wide range study of brain connections.

In this scenario, the multimodal integration of evidences provided by different techniques is the most reliable approach for the investigation of the human connectome. Confirming the reliability of the different approaches and integrating structural and functional multimodal evidences, from both in- and ex-vivo studies, provide more reliable data as well as some technical and conceptual refinements of each approach. This virtuous loop already allowed:

#### REFERENCES


Each of the 9 articles presented in this Research Topic is fully in line with fascinating current debates on the definition, organization, terminology, and conceptualization of the anatomy of human WM.

#### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.


**Conflict of Interest:** 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.

Copyright © 2019 Sarubbo and Petit. 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.

# A Connectomic Analysis of the Human Basal Ganglia Network

Alberto Cacciola<sup>1</sup> \* † , Alessandro Calamuneri 2†, Demetrio Milardi 1, 2, Enricomaria Mormina<sup>2</sup> , Gaetana Chillemi <sup>2</sup> , Silvia Marino<sup>1</sup> , Antonino Naro<sup>1</sup> , Giuseppina Rizzo<sup>2</sup> , Giuseppe Anastasi <sup>2</sup> and Angelo Quartarone1, 2 \*

*1 IRCCS Centro Neurolesi "Bonino Pulejo", Messina, Italy, <sup>2</sup> Department of Biomedical, Dental Sciences and Morphological and Functional Images, University of Messina, Messina, Italy*

The current model of basal ganglia circuits has been introduced almost two decades ago and has settled the basis for our understanding of basal ganglia physiology and movement disorders. Although many questions are yet to be answered, several efforts have been recently made to shed new light on basal ganglia function. The traditional concept of "direct" and "indirect" pathways, obtained from axonal tracing studies in non-human primates and post-mortem fiber dissection in the human brain, still retains a remarkable appeal but is somehow obsolete. Therefore, a better comprehension of human structural basal ganglia connectivity *in vivo*, in humans, is of uttermost importance given the involvement of these deep brain structures in many motor and non-motor functions as well as in the pathophysiology of several movement disorders. By using diffusion magnetic resonance imaging and tractography, we have recently challenged the traditional model of basal ganglia network by showing the possible existence, in the human brain, of cortico-pallidal, cortico-nigral projections, which could be mono- or polysynaptic, and an extensive subcortical network connecting the cerebellum and basal ganglia. Herein, we aimed at reconstructing the basal ganglia connectome providing a quantitative connectivity analysis of the reconstructed pathways. The present findings reinforce the idea of an intricate, not yet unraveled, network involving the cerebral cortex, basal ganglia, and cerebellum. Our findings may pave the way for a more comprehensive and holistic pathophysiological model of basal ganglia circuits.

Keywords: connectivity, connectome, constrained spherical deconvolution, diffusion magnetic resonance imaging, neural pathways, white matter, tractography

#### INTRODUCTION

The current view of basal ganglia neurophysiology is largely based on data from rodents and nonhuman primates using various tract-tracing methods for the study of monosynaptic connections combined with immunohistochemistry and in situ hybridization (Albin et al., 1989; DeLong, 1990). According to this view, the cerebral cortex is widely connected with the basal ganglia via two major projection systems; the so-called "direct" and "indirect" striatofugal pathways originating from segregated populations of striatal projection neurons that exert opposite effects upon basal ganglia outflow. In this model, cortical efferents terminate on specific subcortical areas, which in turn project to the same cortical area that started the initial impulse. The main evidence of this morpho-functional segregation was originally provided for the dorsal striatum, which was classically considered as the main input station in the basal ganglia circuitry (Albin et al., 1989).

#### Edited by:

*Silvio Sarubbo, Ospedale Santa Chiara, Italy*

#### Reviewed by:

*Marina Bentivoglio, University of Verona, Italy Michela Ferrucci, University of Pisa, Italy*

#### \*Correspondence:

*Alberto Cacciola alberto.cacciola0@gmail.com Angelo Quartarone aquartar65@gmail.com*

*† Shared first authorship.*

Received: *08 May 2017* Accepted: *11 September 2017* Published: *26 September 2017*

#### Citation:

*Cacciola A, Calamuneri A, Milardi D, Mormina E, Chillemi G, Marino S, Naro A, Rizzo G, Anastasi G and Quartarone A (2017) A Connectomic Analysis of the Human Basal Ganglia Network. Front. Neuroanat. 11:85. doi: 10.3389/fnana.2017.00085*

**7**

In addition to the striatum, a glutamatergic hyper-direct pathway connecting the cerebral cortex and subthalamic nucleus (STN) has been described in the monkey (Nambu et al., 2002) and would be responsible for a rapid signal transmission from the cortex to the internal segment of the globus pallidus (GPi) (Nambu et al., 2000).

Classically, the anatomical description of the basal ganglia circuitry has been based on axonal tracing in non-human primates and post-mortem fiber dissection in the human brain. However, in the last 20 years, Diffusion Tensor Imaging (DTI) tractography has been used as an anatomical "virtual dissector" for tracking in vivo and non-invasively neural connectivity in the human brain (Basser et al., 1994; Henderson, 2012; Nunnari et al., 2014), by measuring anisotropy diffusion of water molecules (Mori and Van Zijl, 2002).

A DTI study in humans demonstrated that projections from pre-supplementary motor area (pre-SMA) reach only the anterior and middle part of the caudate nucleus and putamen, while the SMA and primary motor area (M1) send dense efferents to the posterior part of putamen and a few to its middle part (Lehéricy et al., 2004a). Further investigations demonstrated that, in humans, the caudate head is connected with several frontal areas (Lehéricy et al., 2004b; Draganski et al., 2008). On the other hand, the caudate tail is technically difficult to study in humans because of its poor detectability with neuroimaging techniques. Studies on primates suggest a major connectivity with temporal areas (Yeterian and Van Hoesen, 1978). This topographical distribution of inputs appears to be maintained in striatal efferents in macaques (François et al., 1994) and in pallidal and nigral efferents in monkeys (Cebus Apella) (Hoover and Strick, 1999).

The existence of the hyper-direct pathway in humans was first demonstrated in a study assessing the STN connectivity with the motor, associative, and limbic brain areas, based on structural and functional connectivity analysis (Brunenberg et al., 2012). The STN is connected with the primary motor cortex, premotor and supplementary motor cortex, and the somatosensory cortex, in lines with studies in non-human primates (Nambu et al., 2002; Brunenberg et al., 2012; Milardi et al., 2015).

In a recent study, we have challenged the traditional model of the basal ganglia network by showing evidence of cortico-pallidal projections in humans by means of Constrained Spherical Deconvolution (CSD) based tractography, which allows to successfully resolve multiple fiber populations insisting over the same voxel and to overcome most of the DTI model issues.

Indeed, we identified a direct cortico-pallidal pathway connecting both the GPi and the external segment of the globus pallidus (GPe) with several cortical regions and running through the internal capsule (Milardi et al., 2015). This pathway may represent an additional route for upstream cortical regulation on basal ganglia circuitry, bypassing the striatum and paralleling the hyperdirect pathway (Smith and Wichmann, 2015).

Furthermore, we showed tractographic evidence of a complementary direct cortico-substantia nigra (SN) pathway paralleling the direct, indirect, hyperdirect, and cortico-pallidal systems (Cacciola et al., 2016a). Several studies revealed corticonigral connections in rodents (Sesack et al., 1989), primates (Frankle et al., 2006) and humans (Kwon and Jang, 2014; Cacciola et al., 2016a). Although in animals it was shown that this direct pathway may represent a glutamatergic cortical regulator of SN activity (Kornhuber et al., 1984), in humans its function and physiological significance are still unknown.

It is worthy to note that diffusion-based tractography is not sufficient to provide anatomical evidence of the existence of a specific pathway, if used alone (Jbabdi and Johansen-Berg, 2011). Indeed, it calculates the highest mathematical probability that water diffuses in a given direction, thus inferring the preferential water diffusivity directionality along white matter bundles. In addition, it should be considered that tractographic reconstruction is not able to clearly distinguish mono- from polysynaptic connections as well as the directionality (afferent or efferent) of the signal transmission (Chung et al., 2011; Parker et al., 2013). Nevertheless, tractography is the only method available to study non-invasively anatomical connectivity in the human brain in vivo.

In the present study, we used a robust diffusion signal model, namely the CSD, for tracking and assessing the connectivity patterns and topographical organization of the basal ganglia network.

#### MATERIALS AND METHODS

#### Participants

We recruited 15 human subjects (9 women, mean age 29 years; age range 25–32 years) with no previous history of disease. The study followed the tenets of the Declaration of Helsinki; written informed consent was signed from all included subjects after explanation of the nature and possible consequences of the procedure. The study was approved by the institutional review board of IRCCS Bonino Pulejo, Messina, Italy (Scientific Institute for Research, Hospitalization and Health Care).

#### Data Acquisition

The study was performed with a 3T Achieva Philips scanner equipped with a 32-channels SENSE head coil. In each subject, a structural 3D high-resolution T1 weighted Fast Field Echo (FFE) sequence was acquired using the following parameters: repetition time 25 ms; echo time 4.6 ms; flip angle 30◦ ; FOV 240 × 240 mm<sup>2</sup> ; isotropic voxel size 1 mm. The acquisition time was 6 min. Furthermore, a 3-D high resolution T2 weighted Turbo Spin Echo (TSE) sequence was obtained using the following parameters: repetition time 2,500 ms; echo time 380 ms; FOV 250 × 250 mm<sup>2</sup> ; in plane reconstruction matrix 312 × 312; slice thickness 0.8 mm. The acquisition time was 9 min and 38 s.

The use of 3D TSE sequence allowed to obtain high-resolution images with a relative short acquisition time, as well as to provide a fine representation of the iron loaded nuclei due to T2<sup>∗</sup> effect linked with the use of a very long echo-time.

Furthermore, a Diffusion weighted dataset (DWI) was obtained with a single-shot EPI sequence using the following parameters: repetition time 11,884 ms, echo-time 54 ms, FOV 240 × 240 mm<sup>2</sup> , isotropic voxel size of 2 mm. Thirty-two diffusion encoded volumes were acquired using a b-value of 1,000 s/mm<sup>2</sup> ; in addition, two un-weighted b0 volumes were acquired, one of which with inverted phase-encoding direction for post-acquisition distortion corrections.

#### DWI Pre-processing and Co-registration

All diffusion images were corrected for motion as well as for Eddy Currents distortion artifacts using a combination of topup and eddy FMRIB Software Library (FSL) tools (http://fsl.fmrib.ox. ac.uk/fsl/fslwiki/); rotational part of transformations were later applied to individual gradient directions. High Resolution T1 and T2 images were then co-registered to preprocessed DWIs using a pipeline outlined in Besson et al. (2014): for each subject, cerebrospinal fluid (CSF) probability maps were estimated from b0 image as well as from T1w volume using New Segment SPM8 utility (http://www.fil.ion.ucl.ac.uk/spm/software/spm8/). Then, CSF of structural scan was warped to match CSF estimated from b0 image using FLIRT (FMRIB's Linear Image Registration Tool) and FNIRT (FMRIB's non-linear registration tool). Subsequently, estimated deformation field was applied to T1w volume. In a separate stage, an affine transformation, using FLIRT FSL command, was performed in order to co-register the T2w volume to T1w scan; eventually, the same warping field previously obtained was applied to this map. In this way, we obtained an overlap between structural images and DWIs as accurate as possible.

### Nuclei Segmentation and Cortical Parcellation

SN and STN segmentation was manually performed on T1w and T2w co-registered volumes by a skilled neuroradiologist, as outlined in a previous our work (Milardi et al., 2016a).

Cortical and subcortical reconstruction and volumetric segmentation were performed on co-registered T1 images with the Freesurfer image analysis suite, which is documented and freely available for download online (http://surfer.nmr.mgh. harvard.edu/). Briefly, this processing involves motion correction and averaging (Reuter et al., 2010) of T1 weighted images, removal of non-brain tissue using a hybrid watershed/surface deformation procedure (Ségonne et al., 2004), segmentation of the subcortical white matter and deep gray matter volumetric structures (Fischl et al., 2004), tessellation of the gray matter white matter boundary, automated topology correction (Ségonne et al., 2007), and surface deformation following intensity gradients to optimally place the gray/white and gray/CSF borders at the location where the greatest shift in intensity defines the transition to the other tissue class (Fischl and Dale, 2000). Once the cortical models were completed, parcellation of the cerebral cortex into units with respect to gyral and sulcal structure (Desikan et al., 2006) was performed. Subsequently, the obtained parcellations and segmentations of each subject were visually inspected and, if needed, manually edited.

### Tractography

To model diffusion signal, we used a modified High Angular Resolution Diffusion Imaging (HARDI) technique called nonnegative CSD: this technique estimates the fiber Orientation Distribution Function (fODF) directly from deconvolution of DW signal with a reference single fiber response function (Tournier et al., 2007, 2008). fODF estimation and tractography were performed using MRtrix software (http://www.mrtrix.org).

By using CSD to extract local fiber orientations we could overcome partial volume effects associated with DTI. Higher b-values permit to resolve smaller angles among fibers (Alexander and Barker, 2005; Tournier et al., 2007); on the other hand, they require longer acquisition times thus increasing the probability to find motion and eddy currents related artifacts. Thus, we preferred a lower bvalue in order to obtain a reasonable quality speed tradeoff.

In our study, spherical harmonic degree was fixed equal to six in order to obtain robustness to noise. During tractographic reconstruction, tracking was stopped in one of the following conditions: step size = 0.2 mm, maximum angle = 10◦ , minimal fODF amplitude = 0.15. The latter parameter allowed to obtain more accurate reconstructions avoiding streamlines to enter GM in deep or passing through CSF; indeed, in those regions, estimated fODF amplitudes are lower than such cut-off. This is a more conservative choice with respect to usual standards, since we preferred to underestimate fiber bundles in order to have more consistent reconstructions (Descoteaux et al., 2009; Tournier et al., 2011; Milardi et al., 2017).

Whole brain probabilistic tractography was run by generating 10 million tracts using the above-mentioned criteria for streamline creation and end. At this stage, a moderately dilated version of WM mask was used both as seed and termination mask, to further reduce the risk of implausible streamlines to be generated.

### Connectogram Construction

To create connectograms of right and left subcortical nuclei (caudate, putamen, GP, STN, and SN), the following pipeline was adopted: for each nucleus, MRtrix tckedit command (include option) was employed to filter the tracks of interest. At the same time, each structure obtained from Freesurfer segmentation was extracted and included within the same command (include option) to isolate the specific tracks linking the basal ganglia to the related area. The same process was repeated for all ROIs obtained from Freesurfer. Subsequently, in-house scripts built in Matlab (http://www.mathworks.com/), release 2013, were run to create the connectograms for each nucleus.

It is worth to mention that, for an accurate extraction of streamlines of interest, and to avoid erroneous track assignation to a given structure of the basal ganglia, appropriate region of avoidances (exclude option of tckedit command) were used (Verstynen et al., 2011).

Connectivity density of the pathways of interest has been estimated by the region-to-region number of streamlines. To summarize the distribution of the connectivity density for each reconstructed pathway we calculated the median density (δ) and standard deviation (SD) from individual subject profiles. Analysis of coefficient of variation (COV), defined as the ratio of the SD to the δ (COV = SD/δ), was also carried out to assess inter-subjects variability.

## RESULTS

By using probabilistic CSD tractography, we reconstructed the undirected connectivity patterns for each nucleus of interest with several supra- and sub-tentorial brain regions. The choice of a proper cutoff for networks is matter of debate in the literature (van Wijk et al., 2010) and probabilistic tractography may yield to spurious results (Rubinov and Sporns, 2010). For these reasons, only connections whose δ was on average above 1% have been considered.

Connectograms estimated for caudate nucleus (**Figure 1A**), putamen (**Figure 1B**), GP (**Figure 2A**), SN (**Figure 2B**), and STN (**Figure 3**) are shown for visualization purposes.

In particular, we found that the neostriatum showed an extensive fronto-parietal network including many motor and non-motor regions such as the precentral, postcentral, paracentral, inferior, middle, and superior frontal gyrus as well as the insula, lateral, and medial orbitofrontal cortices (**Table 1**).

Topographically mapping of GP revealed connectivity patterns with sensory-motor regions, the cerebellum and temporal structures such as the hippocampus and the amygdala. In addition, the nigro-pallidal and pallidal-subthalamic pathways emerged from the connectivity analysis (**Table 2**).

We also evaluated the connectivity profile of the SN which showed the strongest connectivity pattern with both ipsi- and contralateral cerebellum, thus confirming the presence of several parallel cerebellar-basal ganglia circuits and reinforcing the hypothesis of a cerebello-nigral pathway in humans (Milardi et al., 2016a). In addition, SN connectivity was directed toward several fronto-temporo-parietal areas such as the precentral, postcentral, paracentral, superior, inferior parietal gyri, hippocampus, and amygdala (**Table 3**).

Finally, STN showed the strongest connections with the SN, probably due to their adjacent anatomical location, as well as many connectivity patterns with fronto-parietal regions involved in the hyperdirect cortico-subthalamic pathway, such as the precentral, postcentral, superior frontal gyri. In addition, in line with previous findings in animals and humans, STN connectivity analysis revealed a strong interaction with the cerebellar cortex (**Table 3**).

We investigated the consistency of density percentages estimated from our subjects by looking at the COV. While the individual results are reported in **Tables 1**–**3**, summary gathered from all structures are shown for a better overview in **Table 4**. It emerged that the most consistent results were obtained for the putamen, both for the left and right based connectomes (**Table 4**). Higher variability for the left GP results was observed with respect to the right one, with an average COV of 0.92 vs. 0.72 respectively. An inverse pattern was instead

FIGURE 3 | Connectograms showing the connectivity density profile of right and left subthalamic nuclei.

observed for SN, with a slight increment of the variability of the right based connectomes over the left ones (0.90 vs. 0.79 on average). The highest variability between subjects was instead observed for STN connectomes, especially in the left side (**Table 4**).

### DISCUSSION

In the present study we investigated the topography and organization of structural connections of the basal ganglia in humans, showing that the cortico-basal ganglia network consists of several, parallel, and segregated pathways. We have also provided further tractographic evidences of the presence of an extensive anatomical network connecting the cerebellum with the basal ganglia, in line with previous studies in animals and humans.

It should be noted that, while tractography is compelling in being applicable in vivo non-invasively in humans, in practice it suffers from well-documented technical limitations. First of all, this technique cannot disentangle monosynaptic from polysynaptic connections, being unable to reveal the existence of synapses or gap junctions, in addition to the inability to determine the directionality (afferent-efferent) of fiber tracts (Chung et al., 2011; Parker et al., 2013). Moreover, it is worthy to note that tractography cannot prove that the reconstructed pathways are anatomically accurate, since it is not able to distinguish discrete axonal pathways, whose diameter is typically less than 10 µm. Hence, tractographic findings should be carefully considered (Jbabdi and Johansen-Berg, 2011).

From the technical point of view, the presence of intravoxel geometric heterogeneity, which is typical of more than 90% of white matter voxels (Jeurissen et al., 2013), has been faced by using CSD diffusion model, which is known to better elucidate fiber pathways than DTI (Tournier et al., 2007).

Although our results are plausible if compared with previous anatomical descriptions of these circuits, tractographic results should be interpreted with caution. Especially when making use of probabilistic algorithms for streamline generation, false positives may be reconstructed. In this regard, to make our tractographic findings more consistent, we decided to increase the cutoff for tracking generation and termination, therefore using a conservative approach (Descoteaux et al., 2009; Milardi et al., 2016a,b; Cacciola et al., 2017a). This was accomplished at the cost of an underestimation issue.

In addition, the tracking of parallel segregated pathways throughout the cortico-subcortical circuits may be limited by the overlapping nature of the basal ganglia. We tried to limit such issue by combination of inclusion ROIs and regions of avoidance (ROAs).

On the other hand it is worthy to note that it has been demonstrated that this technique permits to obtain good macroscopic neuroanatomical information on white matter fiber bundles by reconstructing streamlines structures containing bundles of axons running along the same direction (Mori and Van Zijl, 2002).

TABLE 1 | Summary of connectomic analysis (in percentages) of caudate nucleus and putamen.


*Standard deviation (SD) and coefficient of variation (COV) are reported for each structure.*

In conclusion, tractography is an anatomical technique by which functional significance can be hypothesized. However, it may provide an interesting perspective for studying altered connectivity patterns in neuropsychiatric and movement disorders related to the cortico-basal ganglia-cerebellar connectome.

#### Cortico-Striatal Connectivity

Cortico-striatal pathway has been traditionally described as a rich set of connections between the putamen as well as the caudate head and body with frontal, parietal and, rarely, occipital regions. Our results are substantially in line with previous findings showing that caudate head and body are principally connected with frontal areas such as the ventral prefrontal cortex, superior frontal gyrus (SFG), and rostral middle frontal gyrus (rMFG) (Lehéricy et al., 2004b). In addition, we here showed orbitofrontal-striatal connectivity, but, in agreement with previous findings (Draganski et al., 2008), no suprathreshold connections were found for the paracentral lobule and premotor area with caudate head and body.

On the other hand, the present findings of a strong connectivity between the precentral gyrus and putamen are in



*Standard deviation (SD) and coefficient of variation (COV) are reported for each structure.*

line with previous data showing that connections with motor areas preferentially involve the putamen rather than the caudate nucleus (Lehéricy et al., 2004a).

These results suggest that the head and body of caudate nucleus are mainly involved in the so-called "associativecognitive loop," being preferentially connected with areas involved in self-awareness (SFG) (Goldberg et al., 2006), in cognitive processes (dlPFC) (Cieslik et al., 2013), in memory retrieval and language (pars opercularis, pars orbitalis, and pars triangularis) (Kostopoulos and Petrides, 2003; Friederici, 2011).

In addition, a functional striatal connectivity with several frontal areas has been demonstrated, suggesting that the dorsal caudate nucleus is more closely related to the lateral orbitofrontal cortex (OFC) than to medial portion (Di Martino et al., 2008). In line with these findings, the structural connectivity analysis of the present study indeed revealed stronger caudate connectivity patterns with the lateral OFC than with the medial one.

On the other hand, we here confirmed that the putamen receives its major input from parietal areas linked with sensorimotor functions. In addition, putamen is connected with some different associative regions such as the superior temporal gyrus (STG), supramarginal and lingual gyri, and insula. In particular, STG connections with putamen may be considered substantially as the human analog of those described in Cebus Apella monkeys using HSV strains as retrograde fiber tracer (Middleton and Strick, 1996). Connections between putamen and pericalcarine and lingual gyri and between putamen and left superior parietal lobule have been described in human brain in a combined resting state-functional MRI and DTI study, showing that FA value changes in these pathways correlated with cognitive processing speed in healthy humans (Ystad et al., 2011).

Although similar considerations can be made for the functional connectivity of putamen with sensorimotor cortices and with insula, we were not able to match the significant functional connectivity reported for anterior cingulate cortex with both caudate and putamen (δ < 1%) (Di Martino et al., 2008).

#### Cortico-STN Connectivity

Although the STN is considered one of the relevant nodes of the "indirect" pathway, it also receives a direct input from the cerebral cortex (Kitai and Deniau, 1981; Nambu et al., 2000). Indeed, Nambu et al. (2002) described a glutamatergic direct cortical input on the STN, conveying excitatory stimuli from motor, associative and limbic brain areas on the GP, bypassing the indirect inhibitor circuit. In line with these findings, by using 7T MRI and tractography, it has been demonstrated that the posterior medial frontal cortex-STN white matter tract strength predicts inter-individual efficacy in stopping a movement in a motor no-go task (Forstmann et al., 2012). Herein, we showed that the connectivity density of the STN is mainly distributed to the precentral, postcentral gyri, and the paracentral lobule, which are cortical areas involved in sensorimotor functions. Moreover, the connectivity pattern linking STN and SFG detected in the present study is in line with the recent idea of a "cognitive STN" involved in decision-making (Weintraub and Zaghloul, 2013).

TABLE 3 | Summary of connectomic analysis (in percentages) of substantia nigra and subthalamic nucleus.


*Standard deviation (SD) and coefficient of variation (COV) are reported for each structure.*

### Cortico-Pallidal Connectivity

Although the cortical inputs are traditionally thought to reach the GP via the direct, indirect, and hyperdirect pathways, there is growing evidence, in animal studies, of the possible existence of a direct route connecting the GP and the cerebral cortex reciprocally. In an anterograde tract tracing study, Naito and Kita (1994) described the presence of a cortico-pallidal projection (which accounted for about 10% of the corticostriatal pathway density) in rodents, linking the medial and lateral precentral cortices to the GPe. More recently, it has been reported that in turn cholinergic and GABAergic neurons within the GPe send direct efferent to the cortex (Chen et al., 2015; Saunders et al., 2015). This suggests a possible major revision of basal ganglia circuits in which GP projections toward the frontal cortex may modulate the activity of different cortical areas. Furthermore, Milardi et al. (2015) highlighted the possible existence of a direct cortico-pallidal pathway in humans based on CSD-tractography; its function, however, remains to be clarified.

In the present study, we provide, for the first time, a quantitative characterization of the GP structural connectivity with the cerebral cortex. Our connectivity analysis showed an anatomical pallidal-temporal network involving the hippocampus and amygdala**,** as well as a sensorimotor-pallidal network involving mainly the precentral, postcentral, and paracentral gyrus, in addition to high-order functions-related areas such as the SFG.

Taking into account that our approach used the caudate, putamen, and thalamic nuclei as exclusion masks in order to avoid the reconstruction of the short subcortical connections within the subcortical basal ganglia network, the described connectivity patterns might reflect the direct information flow between the cerebral cortex and GP, in line with previous animal and human studies (Naito and Kita, 1994; Chen et al., 2015; Milardi et al., 2015; Saunders et al., 2015).

In addition, considering the combination of ROIs and ROAs employed for tract reconstruction, it is likely that the corticopallidal pathway constitutes an additional system separate from the cortico-spinal tract and cortico-pontine pathway. Indeed, Smith and Wichmann (2015) in an editorial note suggested the cortico-pallidal route as an additional regulatory



glutamatergic pathway in basal ganglia circuitry in rodents, monkeys, and humans.

#### Cortico-Nigral Connectivity

Only a few studies have investigated so far the structural connectivity of SN in the human brain using DTI. Two nigral regions have been identified on the basis of connectivity patterns (Menke et al., 2010): an internal region **(**likely corresponding to the SNc) which was mainly connected with the striatum, GP, anterior thalamus, and prefrontal cortex, as well as an external region showing major connectivity with posterior thalamus, ventral thalamus, and motor cortex, which likely corresponded to the SNr, in agreement with the anatomical knowledge (Parent and Hazrati, 1994). More recently Kwon and Jang (2014) evaluated differences in structural connectivity between SN and ventral tegmental area, showing that SN has an extensive network with many frontal, parietal, temporal areas as well as with the cerebellum. We have recently hypothesized the possible existence of a direct cortico-nigral pathway connecting SN with many cortical areas related to motor functions, such as precentral gyrus, postcentral gyrus, and paracentral lobule as well as to the SFG in the prefrontal cortex (Cacciola et al., 2016a).

In the present study the whole brain-based connectivity analysis revealed that the SN is mainly connected with the precentral, postcentral gyri, and paracentral lobules, in addition to frontal and parietal areas involved both in sensorimotor and high order functions. If, on the one hand, the connectivity patterns between prefrontal cortex and SN are consistent with previous findings in primates (Frankle et al., 2006), the sparse connectivity described between OFC and SN in primates was not found to be sufficiently dense in the present study, in line with our previous results (Cacciola et al., 2016a). The presence of an intricate network between the cerebral cortex and SN might play a very important role in the basal ganglia circuitry, as it may play as an upstream control on basal ganglia response by modulating the nigro-striatal system.

It is worthy to note that beyond the strong connectivity between SN and STN, which is likely to be heavily influenced by their close anatomical location, SN revealed consistent connections with the cerebellum, which are the subject of the next section.

#### Cerebellum-Basal Ganglia Connectivity

Last but not least, we here highlighted extensive parallel subcortical loops between the basal ganglia and the cerebellum. Several studies have suggested that cerebellum and basal ganglia are strongly interconnected in physiological and pathological conditions.

In line with previous studies (Bostan et al., 2010, 2013), we here provided further support to the existence of a pathway running between the STN and cerebellar cortex. Indeed, we found that STN showed a strong connectivity profile with the ipsilateral cerebellar cortex. The weaker connectivity observed for the contralateral connections might be due the longer distance covered by these pathways to reach their final targets. Although from tractographic reconstruction we cannot infer about directionality, this pathway may correspond to STN efferents toward the cerebellar cortex demonstrated in monkeys using retrograde virus-tracing technique (Bostan et al., 2010).

In the present study, the whole brain-based connectivity analysis revealed that, at the level of the basal ganglia, the cerebellum interacts not only with the STN but with the GP as well.

Our results show dense connectivity between SN and both ipsilateral and contralateral cerebellar cortex. Our findings are in line with previous fMRI study which demonstrated a bilateral connectivity between SNc and cerebellar vermis (Tomasi and Volkow, 2014) as well as posterior and lateral cerebellar cortex (Zhang et al., 2016). Indirect evidence for a structural-functional interaction between SN and cerebellum comes also from studies in animals which suggested that the focal electrical stimulation of the dentate and fastigial nuclei leads to fluctuation in dopamine turnover in the SN (Nieoullon et al., 1978).

#### CONCLUDING REMARKS

The last decade has been characterized by the growing idea that, in addition to the direct, indirect, and hyperdirect pathways, several other circuits can contribute to modulate the basal ganglia functioning. Herein, we performed a connectivity analysis on diffusion MRI data in order to characterize the topographical distribution of the connections involved in the cortico-basal ganglia-cerebellar network.

Besides confirming the dense striatal and subthalamic structural connectivity with several frontal and parietal areas, we here showed that the GP and SN are densely interconnected with the cerebral cortex as well. Although we cannot provide any direct evidence on the functional significance of these connections, the findings lead to speculate that the sensorimotorpallidal pathway plays a relevant role in motor control within the cortico-basal ganglia-cortical loop, allowing for a faster information flow between the cerebral cortex and the basal ganglia with respect to the direct, indirect and hyperdirect pathways, as previously hypothesized (Cacciola et al., 2016b).

In addition, our analysis revealed that the STN, GP, and SN are structurally interconnected with the cerebellum. These findings reinforce the idea that the basal ganglia are part of an extensive intricate network with the cerebellum thus mediating the afferent and efferent information flow to integrate inputs of different modalities and to rapidly provide an adaptive motor response (Cacciola et al., 2017b).

Finally, it is worth to note that the interplay between the basal ganglia and cerebellum may be relevant in several movement disorders such as dystonia and Parkinson's disease considering the emergent idea of a compensatory and modulatory function of the cerebellum either at the early stages of the disease or during its progression (Wu and Hallett, 2013; Neumann et al., 2015; Chillemi et al., 2017).

#### REFERENCES


Our findings may pave the way for a more comprehensive and holistic pathophysiological model of basal ganglia circuits.

#### AUTHOR CONTRIBUTIONS

AlbC: Study concepts/study design, data analysis, data interpretation, literature research, and draft the manuscript. AleC: Study concepts/study design, data analysis, and draft the manuscript. DM: Data acquisition, data interpretation, and literature research. EM, GC, and AN: Data analysis, data interpretation, and literature research. SM: Data acquisition and data interpretation. GR: Data interpretation and literature research. GA: Guarantor of integrity of entire study, data interpretation, and manuscript revision for important intellectual content. AQ: Study concepts/study design, guarantor of integrity of entire study, and manuscript revision for important intellectual content. All the authors approved the final version of the manuscript.


interindividual efficacy in stopping a motor response. Neuroimage 60, 370–375. doi: 10.1016/j.neuroimage.2011.12.044


**Conflict of Interest Statement:** 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.

Copyright © 2017 Cacciola, Calamuneri, Milardi, Mormina, Chillemi, Marino, Naro, Rizzo, Anastasi and Quartarone. 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) or licensor 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.

## The Superior Fronto-Occipital Fasciculus in the Human Brain Revealed by Diffusion Spectrum Imaging Tractography: An Anatomical Reality or a Methodological Artifact?

#### Yue Bao, Yong Wang, Wei Wang and Yibao Wang\*

Department of Neurosurgery, The First Affiliated Hospital of China Medical University, Shenyang, China

The existence of the superior fronto-occipital fasciculus (SFOF) in the human brain remains controversial. The aim of the present study was to clarify the existence, course, and terminations of the SFOF. High angular diffusion spectrum imaging (DSI) analysis was performed on six healthy adults and on a template of 842 subjects from the Human Connectome Project. To verify tractography results, we performed fiber microdissections of four post-mortem human brains. Based on DSI tractography, we reconstructed the SFOF in the subjects and the template from the Human Connectome Project that originated from the rostral and medial parts of the superior and middle frontal gyri. By tractography, we found that the fibers formed a compact fascicle at the level of the anterior horn of the lateral ventricle coursing above the head of caudate nucleus, medial to the corona radiate and under the corpus callosum (CC), and terminated at the parietal region via the lower part of the caudate nucleus. We consider that this fiber bundle observed by tractography is the SFOF, although it terminates mainly at the parietal region, rather than occipital lobe. By contrast, we were unable to identify a fiber bundle corresponding to the SFOF in our fiber dissection study. Although we did not provide definite evidence of the SFOF in the human brain, these findings may be useful for future studies in this field.

#### Edited by:

Silvio Sarubbo, Ospedale Santa Chiara, Italy

#### Reviewed by:

Basilis Zikopoulos, Boston University, United States Ilyess Zemmoura, François Rabelais University, France Emmanuel Mandonnet, Lariboisière Hospital, France

#### \*Correspondence:

Yibao Wang cmuwyb@hotmail.com

Received: 17 April 2017 Accepted: 27 November 2017 Published: 13 December 2017

#### Citation:

Bao Y, Wang Y, Wang W and Wang Y (2017) The Superior Fronto-Occipital Fasciculus in the Human Brain Revealed by Diffusion Spectrum Imaging Tractography: An Anatomical Reality or a Methodological Artifact? Front. Neuroanat. 11:119. doi: 10.3389/fnana.2017.00119 Keywords: superior fronto-occipital fasciculus, diffusion spectrum imaging, tractography, white matter, fiber dissection

### INTRODUCTION

The superior fronto-occipital fasciculus (SFOF) was first described as a longitudinal ridge in autopsies of acallosal patients by Eichler (1878). Since then, numerous post-mortem studies have proposed a range of descriptions, hypotheses and terminology of the SFOF. Sachs (1892) made the important observation that the SFOF of acallosal patients was composed of callosal fibers that failed to pass through the interhemispheric midline, which was verified in later more advanced studies (Forkel et al., 2014). As well as confirming that the SFOF represented an association bundle connecting frontal and occipital regions in the acallosal brain, Onufrowicz (1887) and Dejerine (1895) extended these findings to the healthy brain. Thus, the prevailing consensus proposed by Dejerine is that the SFOF is a distinct 2–3 mm wide white matter located medial to the corona radiate, superolateral to the caudate nucleus, and ventrolateral to the corpus callosum (CC), and which interconnects the frontal and occipital lobes as association fibers.

Using isotope tract-tracing in the monkey brain, Schmahmann and Pandya confirmed that the SFOF courses above the caudate nucleus, medial to the corona radiate, and lateral and ventral to the CC, connects the parietal lobe with the frontal lobe, and plays a role in the rapid top-down modulation of visual processing and spatial aspects of cognitive processing (Bar et al., 2006; Schmahmann and Pandya, 2007). Axonal tracing results of non-human primates are generally considered the ''gold standard'' guide for human white matter connective patterns (Thiebaut de Schotten et al., 2012). However, with evolution, the frontal and parieto-occipital regions of the human brain differ markedly from their monkey counterparts. Thus, extrapolating results from monkey studies to humans can be misleading (Rilling et al., 2008). Indeed, despite being well described in human, the inferior fronto-occipital fasciculus (IFOF) does not exist in the monkey brain (Schmahmann et al., 2007; Forkel et al., 2014).

In the last decade, the development of diffusion tensor imaging (DTI) and diffusion tractography have provided a unique opportunity to study white matter architecture in vivo (Dick and Tremblay, 2012). Indeed, a number of fiber tractography DTI studies have identified the SFOF in the human brain (Catani et al., 2002; Wakana et al., 2004; Makris et al., 2007). By contrast, using tractography with DTI Meola et al. (2015) did not identify any tracts in the previously described SFOF location. Further, the SFOF was not identified by tractography with a diffusion imaging approach based on spherical deconvolution (Forkel et al., 2014), and the ''SFOF'' was proposed to represent a branch of the superior longitudinal fasciculus (SLF). Recent microsurgical anatomical studies in post-mortem human brains also reported that the SFOF is actually the superior thalamic peduncle (STP), rather than association fibers (Türe et al., 1997; Meola et al., 2015).

These contrasting fiber dissection and in vivo tractography findings may be attributed to the methodological limitations of these techniques. For example, DTI-based tractography is unable to resolve multiple fiber orientations within an magnetic resonance imaging (MRI) voxel, and as such, cannot resolve fiber crossings either as tract intersections in the white matter or the intricate architecture of the gray matter (Mori and van Zijl, 2002), resulting in artifacts and false tracts (Le Bihan et al., 2006; Fernandez-Miranda, 2013; Wang et al., 2013). In recent years, high-angular-resolution diffusion imaging (HARDI) has been used to overcome the deficiencies of DTI. HARDI is based on diffusion spectrum imaging (DSI), and is reconstructed by generalized Q-sampling imaging (Yeh et al., 2010). Thus, in the present study we mapped the SFOF in six healthy subjects, as well as using a template approach (HCP-842 Atlas), using HARDI tractography. We also performed a fiber dissection study with cross-validation to tractography findings. Our findings provide further clarification of the anatomy of the SFOF, and may be beneficial for planning approaches for neurosurgical procedures.

#### MATERIALS AND METHODS

### Participants

Six neurologically healthy adults (three men; all right handed; age range: 26–43 years) were included in this study. All procedures were approved by the Ethics Committee of the First Affiliated Hospital of China Medical University, and written consent was obtained from all participants before testing. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the the Ethics Committee of the First Affiliated Hospital of China Medical University.

We also performed fiber tracking on a template of 842 subjects from the Human Connectome Project (HCP-842), which represents the largest and highest quality data set available to date. The Wu-Minn HCP consortium is an ongoing project led by Washington University, University of Minnesota and Oxford University, involving a systematic effort to map macroscopic human brain circuits and their relationship to behavior in a large population of healthy adults. Data from a total of 900 subjects were released for the first four quarters (Q1–Q4, 2015), and 842 subjects had diffusion scans. The reconstructed data of the 842 subjects were averaged to create a representative template (the HCP 842-subject template is freely downloadable at: http://dsi-studio.labsolver.org/download-images).

#### Image Acquisition and Reconstruction

For the subject-specific approach, DSI data were acquired on 3T Tim Trio System (Siemens) using a 32-channel coil. This involved a 43-min, 257-direction scan using a twice-refocused spin-echo EPI sequence and multiple q values (Wedeen et al., 2005; repetition time [TR] = 9.916 ms, echo time [TE] = 157 ms, voxel size = 2.4 × 2.4 × 2.4 mm, field of view = 231 × 231 mm, bmax = 7000 s/mm−<sup>2</sup> ). We also included the high-resolution anatomical imaging as anatomical comparisons, employing a 9-min T1-weighted axial magnetization prepared rapid gradient echo sequence (TE = 2.63 ms, TR = 2110 ms, 176 slices, flip angle = 8◦ , field of view = 256 × 256 mm−<sup>2</sup> , voxel size = 0.5 × 0.5 × 1.0 mm<sup>3</sup> ). DSI data were reconstructed using a generalized Q-sampling imaging approach (Yeh et al., 2010). The orientation distribution functions were reconstructed to 362 discrete sampling directions and a mean diffusion distance of 1.2 mm.

#### The Human Connectome Project 842-Subject Template

The HCP-842 Atlas was constructed using diffusion MRI data from a total of 842 subjects from the HCP (2015 Q4, 900-subject release). The 842 subjects underwent diffusion scans in a Siemens 3T Skyra scanner using a two dimensional spin-echo single-shot multiband EPI sequence with a multi-band factor of three and a monopolar gradient pulse (Sotiropoulos et al., 2013). The spatial resolution was 1.25 mm isotropic (TR = 5500 ms, TE = 89 ms). A multishell diffusion scheme was used for acquiring the diffusion images, with b-values of 1000, 2000 and 3000 s/mm<sup>2</sup> . The total number of diffusion sampling directions was 270. The total scanning time was approximately 55 min (Van Essen et al., 2013). Diffusion data were collected using a Q-sampling imaging approach.

#### Fiber Tracking and Analysis

DSI Studio (open-source diffusion MRI analysis tool, freely downloadable at http://dsi-studio.labsolver.org) was used for fiber tracking. A whole brain seeding approach using multiple region of interest (ROI) masks was performed. In voxels with multiple fiber orientations, fiber tracking was initiated separately for each orientation, and fiber progression continued with a step size of 1.2 mm, minimum fiber length of 20 mm and turning angle threshold of 60◦ . If multiple fiber orientations existed in the current progression location, the fiber orientation that was nearest to the incoming direction and formed a turning angle smaller than 60◦ was selected to determine the next moving direction (Fernandez-Miranda et al., 2015). To smooth each track, the next moving directional estimate of each voxel was weighted by 20% of the previous incoming direction and 80% of the nearest fiber orientation. This progression was repeated until the quantitative anisotropy of the fiber orientation dropped below a preset threshold (0.03–0.06 depending on the subject) or there was no fiber selected within the 60◦ angular range in the progression (Yeh et al., 2013). Once tracked, all streamlines were saved in the TrackVis file format. For comparison, FreeSurfer<sup>1</sup> was used to segment cortical gyral ROIs based on previous brain atlases using each participant's T1-weighted magnetization prepared rapid axial gradient echo image (Desikan et al., 2006).

According to prior anatomical studies reported by Dejerine (1895) and experimental animal studies (Schmahmann and Pandya, 2007), the SFOF is a distinct fiber bundle coursing between the ependyma of the lateral ventricle and the internal capsule, superolateral to the caudate nucleus, medial to the corona radiate, and under the CC. Using this information, we tried to reconstruct the SFOF. Three ROIs were drawn on the coronal quantitative anisotropy color map to select the anteroposterior direction fibers that pass through all three ROIs. On the sagittal plane, the rostral ROI was located at the level of the anterior commissure, the middle ROI was located at the level of the thalamus, and the caudal ROI was located at the level of the pineal. On the coronal plane, all three ROIs were located superolateral to the caudate nucleus, medial to the internal capsule and under the CC (**Figure 1**).

We also reconstructed the SLF-II, IFOF and the cingulum bundle (CB) to analyze the spatial relationship of the SFOF with these association fibers. For the SLF, the angular gyrus, supramarginal gyrus, superior parietal lobule, and precuneus were used as seed regions, and an ROI mask was drawn along the precentral region on the coronal plane (Wang et al., 2016). For the IFOF, the frontal lobe was selected as the seeding region, and two ROIs were drawn at the level of the central sulcus and ventral part of the external capsule on the coronal plane (Wu et al., 2016a). For the CB, the cingulate cortex was selected as the seeding region (Wu et al., 2016b).

### Fiber Dissection Technique

Eight hemispheres from four normal brains (36–70 years old, three men) were obtained at routine autopsy. Operations were approved by the ethics committee at China Medical University. The fiber dissection study was performed using the Klingler technique (Agrawal et al., 2011). The specimens were fixed in a 10% formalin aqueous solution for at least 4 weeks, and we then removed the vessels, arachnoid and pia-matter, which were frozen for an additional 2 weeks at −16◦C. Because of the freezing progress, water crystallization disrupted the structure of the gray matter, allowing us to peel off the cortex from brain surface and isolate the fiber bundles in their glial sheets. We performed fiber dissection studies at the Surgical Neuroanatomy Lab of China Medical University using microsurgical instrumentation and a surgical microscope (6–40 magnification; OPMI CS-NC; Carl Zeiss).

#### RESULTS

#### Fiber Dissection Findings

We performed a fiber dissection from the lateral to medial direction. After the arcuate fasciculus (AF), SLF-II and putamen were removed, the internal capsule was exposed (**Figure 2A**). We then removed the anterior part of the internal capsule, and the head of caudate and the thalamus were then exposed (**Figure 2B**). The thalamus was medial to the medial and posterior part of the internal capsule, further remove the internal capsule allowing us to visualize the anterior thalamic peduncle (ATP), superior thalamic peduncle (STP) and posterior thalamus peduncle (PTP; **Figure 2C**). After removal of the internal capsule and thalamus, we were able to expose the most medial part of thalamic peduncles (**Figures 2D,E**). After removing the whole thalamus, we observed a few fibers coming from the parietal region, which had a close relationship to the stria terminalis (ST; **Figure 2F**). However, in our fiber dissection, we did not observe any fibers from the frontal lobe that curved around the upper part of the head of caudate and coursed to the parietal and occipital region (as observed by fiber tractography; **Figure 3**).

We also performed a fiber dissection from the medial to lateral direction. After removal of the cortex, the hippocampus, the CC, and the fornix, the caudate and the thalamus were exposed (**Figure 4A**). We then removed the caudate, allowing us to see the internal capsule and ST (**Figure 4B**). Further dissection was performed to identify fibers from the frontal lobe. After we removed the ependyma, the anterior portion of the subcallosal stratum, and the radiation of the CC, we found a few fibers that traversed from the frontal lobe to the superior part of the thalamus, which were likely part of the STP (**Figure 4C**).

During our fiber dissection, we were unable to find a bundle of fibers that coursed around the head of caudate and traveled to the parietal and occipital regions. Above the head of the caudate and medial to the internal capsule, we did observe some

<sup>1</sup>http://surfer.nmr.mgh.harvard.edu

fibers coming from the frontal lobe (**Figure 5A**, arrow). However, further dissection found that the fibers coursed inferior to the thalamus, and were likely part of the STP (**Figure 5B**).

### Trajectory of the SFOF in the Human Brain

#### Subject-Specific Tractography Findings

In the fiber tractography study, we found an anteroposterior bundle of fibers passing through the ROIs in 12 hemispheres (**Figure 6**), which mainly connected the prefrontal cortex with the parieto-occipital region. Based on the cortical parcellation, we performed analysis of the cortical termination of the SFOF (Desikan et al., 2006). This distinct group of fibers mainly originated from the superior frontal gyrus (SFG) and middle frontal gyrus (MFG), and formed a compact fascicle at the level of the anterior horn of the lateral ventricle. The fibers ran posteriorly under the CC, medial to the internal capsule and lateral to the lateral ventricle, and then ran postero-superiorly, as they terminated mainly at the superior parietal lobule and precuneus (**Figure 7**).

#### Template Tractography Findings

The DSI template (HCP-842) showed similar results to the subject-specific findings above (**Figure 8**).

### Spatial Relationship of the SFOF with Adjacent Association Tracts

The IFOF, SLF-II and CB were reconstructed to show the spatial relationship of the SFOF with adjacent association tracts (**Figure 9**). The IFOF passed from the prefrontal region to the occipital region, and the external/extreme capsule level fibers of the IFOF narrowed into a compact fascicle, forming the stem of the IFOF (Wu et al., 2016a). We also simultaneously reconstructed the IFOF and SFOF in one hemisphere, and found that the stem of the SFOF was located above the stem of the IFOF (**Figure 9A**).

The CB is a critical white matter fiber tract in the brain, which forms connections between the frontal lobe, parietal lobe, and the temporal lobe (Wu et al., 2016b). We simultaneously reconstructed the CB and SFOF in one hemisphere, and found that the CB was a sickle-shaped association bundle that nearly encircled the CC and extended to the temporal lobe, and that the CB was located above and medial to the SFOF (**Figure 9B**).

The SLF-II is a component of the SLF that runs on the lateral aspect of the hemisphere and connects the parietal lobe with the caudal MFG and the dorsal precentral gyrus. We simultaneously reconstructed the SLF-II and SFOF in one hemisphere, and

FIGURE 2 | Step by step fiber dissection from lateral to medial of right hemisphere. (A) The arcuate fasciculus (AF), superior longitudinal fasciculus (SLF) II and putamen were removed, the internal capsule was exposed. (B) After removing the anterior part of internal capsule, the head of caudate nucleus, thalamus, anterior thalamus peduncle and superior thalamus peduncle were exposed. (C) The anterior, superior and posterior thalamic peduncle. (D) Further removal of the internal capsule and thalamus, we could exposed the most medial part of thalamus peduncle. (E) Close-up and bottom to up view of the dotted rectangle in (D), it is clear that lots of fiber bundles reach thalamus. (F) After removal of the whole thalamus, some fibers as the white arrow shows coming from the parietal region have a close relationship to stria terminalis (ST). Int. Cap, internal capsule; Cor. Rad, corona radiate; Thal, thalamus; Caud, caudate nucleus; STP, superior thalamus peduncle; ATP, anterior thalamus peduncle; PTP, posterior thalamus peduncle; ST, stria terminalis.

FIGURE 3 | (A) Fiber tractography of the SFOF (subject 3, right hemisphere); lateral view of the left SFOF, we found fibers originated from MFG curved around the upper part of the head of caudate nucleus and coursed to SPL. (B) The fibers related to ST. (C) From (C) to (B), we did not find some fiber bundles coursed above the head of caudate nucleus and connected the parietal region as the fiber tractography (dotted line) showed. But we did find some fibers interconnected the parietal lobe with the bottom part of caudate head via stria ternimalis as the white arrow showed. MFG, middle frontal gyrus; SPL, superior parietal lobe; caud, caudate; ST, stria terminalis.

found that the SLF-II was located superolateral to the SFOF (**Figure 9C**).

### DISCUSSION

Since Dejerine (1895) reported that SFOF exits in the health human brain, SFOF has never been identified in recent fiber dissection studies (Türe et al., 1997; Meola et al., 2015). We did not find a bundle of fibers corresponding to the SFOF in our fiber dissection study either. In the ROI, we found fibers coming from the frontal lobe that coursed above the head of caudate and media to the internal capsule. However, further dissection demonstrated that these fibers belonged to the STP. During our lateral to medial fiber dissection, after removal of the whole thalamus we found some fibers coming from the parietal region that were not the PTP, but which had a close relationship to the ST. However, because of the crossing fibers and complicated structure in this region, we could not identify a distinct bundle of fibers coming from the frontal lobe and connecting with these fibers.

In contrast to our fiber dissection findings, we reconstructed the SFOF in six subjects and in a template of 842 subjects based on the DSI method. We found that the shape and trajectory pattern of the SFOF were similar to previous DTI studies

FIGURE 4 | Fiber dissection from medial to lateral of left hemisphere. (A) Caudate nucleus and thalamus were exposed. (B) After removal of the caudate nucleus, internal capsule and ST were exposed. (C) After removal of the ependyma, anterior portion of the subcallosal stratum and the radiation of the CC, we found the superior thalamus peduncle which came from the frontal lobe and reach the superior part of the thalamus. (D) Fiber tractography of the SFOF (subject 3, right hemisphere). We did not find an affirmative bundle of fibers interconnect the frontal lobe with parietal region as the dotted line in (C) and the SFOF in (D) showed. Thal, thalamus; Caud, caudate nucleus; Int. Cap, internal capsule; ST, stria terminalis; STP, superior thalamus peduncle; MFG, middle frontal gyrus; SPL, superior parietal lobe.

(Catani et al., 2002; Makris et al., 2007). However, recently Meola et al. (2015) reported that prior SFOF reconstructions based on DTI method were generated by false continuations between the STP and ST, and between ST and the PTP. So, can we trust our findings based on DSI tractography? We think we can't fully trust the findings of our DSI tractography study, on the one hand, we did not identify the SFOF in our fiber dissection study; on the other hand, any diffusion tensor tractography may have a false connection due to the biological signals noise, partial volume effects, magnetic resonance resolution and tracking algorithm (Chen and Song, 2008; Mohammadi et al., 2013). However, fiber dissection technique also has limitations for

FIGURE 6 | Diffusion spectrum imaging (DSI) tractography of the SFOF in 12 hemispheres of six subjects on sagittal view. Order: subject 1, 2, 3, 4, 5, 6. Left and right SFOF had a similar location, shape and trajectory in all 12 hemispheres. L, left; R, right.

resolving crossing fibers and following long fibers (Wang et al., 2016). Besides, the individual skills and anatomical knowledge of the investigators may affect the fiber dissection results. We think the results of our DSI tractography study also have some significance. Schmahmann et al. (2007) used DSI and isotope tract-tracing (gold standard) to trace long association fibers in the monkey brain, are found that the two methods showed a remarkable concordance. In our tractography study, to filter out the fibers that belonged to STP and PTP, we set the middle ROI above the thalamus on the sagittal plane, and in the same position as the rostral and caudal ROIs on the coronal plane (**Figure 1C**). Using three ROIs, we reconstructed the SFOF that coursed above the dorsal thalamus (**Figure 4D**), which was distinguishable from the STP-ST-PTP. If the DSI tractographic SFOF really exist, based on cortical parcellation (Desikan et al., 2006), we found that the SFOF would originate from the rostral and medial parts of the superior and MFG (Brodmann area [BA] 8, 9 and 10), the fibers formed a compact fascicle at the level of the anterior horn of the lateral ventricle coursing above the head of caudate, medial to the corona radiate and

FIGURE 9 | Spatial relationship of the SFOF with the adjacent association tracts. (A) Sagittal view showed the SFOF located above the inferior fronto-occipital fasciculus (IFOF). (B) Sagittal view showed the cingulum bundle (CB) was above and medial to the SFOF. (C) Sagittal view showed the SLF-II was above and lateral to the SFOF. SFOF, superior fronto-occipital fasciculus; IFOF, inferior fronto-occipital fasciculus; CB, cingulum bundle; SLF, superior longitudinal fasciculus.

Literally speaking, SFOF is a bundle of fibers connecting frontal and occipital lobes, however, our tractography study and some previous DTI studies (Catani et al., 2002; Makris et al., 2007) show that the SFOF mainly connects frontal lobe with parietal region. If the SFOF really exists in humans, we think that the differences between our result and previous studies of classical neuroanatomy and monkey brain may due to the limitation of gross dissection and brain evolution of human beings.

Despite reconstructing the SFOF based on advanced HARDI in the present tractography study, we cannot conclude that SFOF exists in the human brain as we did not find it by fiber dissection, and there is also potential for false connections in tractography data. The combination of ex vivo DSI scans with fiber dissection in the same human brain, or of in vivo DSI scan with fiber dissection in the same non-human primate brain, the comparison will be more compelling. In addition, new method like polarized light imaging (PLI), which can identify nerve fibers and their spatial orientation in microtome sections of postmortem brains based on the birefringent properties of the myelin sheaths (Dammers et al., 2012; Mollink et al., 2017), can further verify our DSI tractography study.

There were several limitations in our study. First, a subjectspecific approach was completed in only six subjects. More subjects are required to manifest whether the SFOF can be reconstructed in each subject. Second, although fiber dissection can provide key macroscopic information on fiber tract anatomy, it has limited value for the study of white matter trajectory in areas of crossing fibers. Finally, the potential functional role of the SFOF in the present study is only based on correlation analysis, and remains to be proven.

#### CONCLUSION

Using HARDI tractography, we identified the SFOF in six human subjects, and described the trajectory and cortical endpoint of the SFOF. Based on its shape, trajectory, and anatomical relationship with related structures, we consider that this fiber bundle observed by tractography may be the SFOF, although it terminates mainly at the parietal region, rather than occipital lobe. Nevertheless, we were unable to identify a fiber bundle corresponding to the SFOF in our fiber dissection study, which is likely because of limitations in the fiber dissection technique. Although we did not provide definite evidence of the SFOF in the human brain, these findings may be useful for future studies in this field.

#### REFERENCES


#### AUTHOR CONTRIBUTIONS

YiW and YB: conceived and designed the experiments, wrote the article and revised the manuscript. YB, YiW and YoW: performed the experiments. YB and WW: data interpretation and picture preparation. YiW and YoW: collected the biopsy samples. YiW and WW: contributed reagents/materials/analyses tools.

#### FUNDING

This study was supported by grants from the National Natural Science Foundation of China (Grant No. 31540077 to Dr. YiW).


Yeh, F.-C., Wedeen, V. J., and Tseng, W.-Y. I. (2010). Generalized q-sampling imaging. IEEE Trans. Med. Imaging 29, 1626–1635. doi: 10.1109/TMI. 2010.2045126

**Conflict of Interest Statement**: 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.

Copyright © 2017 Bao, Wang, Wang and Wang. 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) or licensor 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.

fnana-11-00138 January 13, 2018 Time: 15:58 # 1

# Evidence for Plastic Processes in Migraine with Aura: A Diffusion Weighted MRI Study

Nikoletta Szabó1,2 \*, Péter Faragó1,2, András Király1,2, Dániel Veréb<sup>1</sup> , Gergo Csete ˝ 1 , Eszter Tóth<sup>1</sup> , Krisztián Kocsis<sup>1</sup> , Bálint Kincses<sup>1</sup> , Bernadett Tuka<sup>3</sup> , Árpád Párdutz<sup>1</sup> , Délia Szok<sup>1</sup> , János Tajti<sup>1</sup> , László Vécsei1,3 and Zsigmond T. Kincses<sup>1</sup> \*

<sup>1</sup> Neuroimaging Research Group, Department of Neurology, Albert Szent-Györgyi Clinical Center, University of Szeged, Szeged, Hungary, <sup>2</sup> Central European Institute of Technology, Brno, Czechia, <sup>3</sup> MTA-SZTE Neuroscience Research Group, Szeged, Hungary

Background: Formerly white matter abnormalities in a mixed group of migraine patients with and without aura were shown. Here, we aimed to explore white matter alterations in a homogeneous group of migraineurs with aura and to delineate possible relationships between white matter changes and clinical variables.

#### Edited by:

Laurent Petit, Centre National de la Recherche Scientifique (CNRS), France

#### Reviewed by:

Filippo Brighina, Università degli Studi di Palermo, Italy Antonio Russo, University of Campania Luigi Vanvitelli, Italy

#### \*Correspondence:

Zsigmond T. Kincses kincses.zsigmond.tamas@med. u-szeged.hu Nikoletta Szabó sznikol@yahoo.com

Received: 10 September 2017 Accepted: 26 December 2017 Published: 17 January 2018

#### Citation:

Szabó N, Faragó P, Király A, Veréb D, Csete G, Tóth E, Kocsis K, Kincses B, Tuka B, Párdutz Á, Szok D, Tajti J, Vécsei L and Kincses ZT (2018) Evidence for Plastic Processes in Migraine with Aura: A Diffusion Weighted MRI Study. Front. Neuroanat. 11:138. doi: 10.3389/fnana.2017.00138 Methods: Eighteen patients with aura, 25 migraine patients without aura and 28 controls were scanned on a 1.5T MRI scanner. Diffusivity parameters of the white matter were estimated and compared between patients' groups and controls using whole-brain tract-based spatial statistics.

Results: Decreased radial diffusivity (p < 0.036) was found bilaterally in the parietooccipital white matter, the corpus callosum, and the cingular white matter of migraine with aura (MwA) patients compared to controls. Migraine without aura (MwoA) patients showed no alteration compared to controls. MwA compared to MwoA showed increased fractional anisotropy (p < 0.048) in the left parieto-occipital white matter. In MwA a negative correlation was found between axial diffusivity and disease duration in the left superior longitudinal fascicle (left parieto-occipital region) and in the left corticospinal tract (p < 0.036) and with the number of the attacks in the right superior longitudinal fascicle (p < 0.048).

Conclusion: We showed for the first time that there are white matter microstructural differences between these two subgroups of migraine and hence it is important to handle the two groups separately in further researches. We propose that degenerative and maladaptive plastic changes coexist in the disease and the diffusion profile is a result of these processes.

#### Keywords: DTI, microstructure, neurodegeneration, plasticity, TBSS, white matter

**Abbreviations:** AD, diffusivity parallel; CSD, cortical spreading depression; DTI, diffusion tensor imaging; FA, fractional anisotropy; FSL, FMRIB Software Library; GLM, general linear model; MD, mean diffusivity; MwA, migraine with aura; MwoA, migraine without aura; RD, perpendicular; TBSS, tract-based spatial statistics.

### INTRODUCTION

fnana-11-00138 January 13, 2018 Time: 15:58 # 2

In 20% of cases, migraine is preceded by focal neurological symptoms, such as visual, sensory, motor or verbal disturbances, called aura. This special form of migraine [MwA (Lauritzen, 1994)] starts with a different trigger than MwoA. CSD, the probable electrophysiological correlate for the aura phenomenon, is thought to occur exclusively in MwA (Lauritzen, 1994; Hadjikhani et al., 2001). Debate continues whether MwA and MwoA are separate entities; some studies question the viability of a distinction between them (Ranson et al., 1991; Manzoni and Torelli, 2008).

Imaging studies revealed that migraine is associated with structural changes affecting gray (Kim et al., 2008; Schmidt-Wilcke et al., 2008; Schmitz et al., 2008; Valfre et al., 2008; May, 2009; Dai et al., 2015) and white matter (Szabo et al., 2012; Chong and Schwedt, 2015). However, there have only been a few studies on how the structural parameters of the two sub-groups of migraine differ from each other (Granziera et al., 2013; Rocca et al., 2014). Most studies investigating white matter integrity in migraineurs either found no significant difference between the two groups, used mixed groups, or assessed just one of them (Rocca et al., 2003, 2006; Chong and Schwedt, 2015; Messina et al., 2015). There are, however, several findings that have brought dissimilarities to light, e.g., the diffusion parameters of the visual pathway (Granziera et al., 2006; Rocca et al., 2008) and the thalamocortical tract (DaSilva et al., 2007). These conflicting reports, the incomplete understanding of mechanisms underlying the disorder and the absence of reliable biomarkers urges further research in the area. Previously we showed the alterations of diffusion MRI measured white matter microstructure in a mixed group of migraine patients with and without aura using TBSS (Szabo et al., 2012). In the current investigation, we extended our former results and went on to assess the white matter microstructural alteration in a group of patients with MwA and MwoA. The main goal of this study was to clarify that microstructural changes of the white matter are distinct within subtypes of migraine and to determine whether these alterations are related to clinical parameters.

### MATERIALS AND METHODS

#### Subjects

Eighteen MwA, 25 MwoA were recruited from outpatients of the Headache Outpatient Clinic of the Department of Neurology. The diagnosis was based on the criteria of the International Headache Society (Headache Classification Committee of the International Headache Society, 2013). Hamilton Depression Scale (Hamilton, 1960) was used to screen for depression and patients with more than eight points were excluded from the study. The patients had no other neurological diseases. Twentyeight controls with no neurological disorders were recruited into the study. The controls were selected to match the age and gender distribution of both patient groups. The Regional Human Biomedical Research Ethics Committee of the University of Szeged approved the study (87/2009). Written informed consent was obtained from all participants.

### MRI Acquisition

All MRI acquisitions were in the interictal period, minimum 1 week after the last headache attack. MRI acquisitions were carried out on a 1.5 T GE Signa Excite HDxt MR Scanner (GE Healthcare, Chalfont St. Giles, United Kingdom). Three dimensional spoiled gradient echo images [FSPGR: echo time (TE): 4.1 ms; repetition time (TR): 10.276 ms; matrix: 256 <sup>∗</sup> 256; field of view (FOV): 25 cm <sup>∗</sup> 25 cm; flip angel: 15◦ ; in-plane resolution: 1 mm <sup>∗</sup> 1 mm; slice thickness: 1 mm] and 60 directions diffusion-weighted images with 6 non-diffusionweighted reference volume [TE: 93.8 ms; TR: 16.000 ms; matrix: 96 <sup>∗</sup> 96; FOV: 23 cm <sup>∗</sup> 23 cm; flip angle: 90◦ ; in-plane resolution: 2.4 mm <sup>∗</sup> 2.4 mm; slice thickness: 2.4 mm; b: 1000 s/mm<sup>2</sup> ; number of excitations (NEX): 2; array spatial sensitivity encoding technique (ASSET) factor = 2] were acquired for all subjects.

### Data Analysis

Image analysis was carried out using tools of FSL (Smith et al., 2004). Diffusion data were corrected for eddy currents and movement artifacts by 12 degrees of freedom affine linear registration to the first non-diffusion-weighted reference image (Jenkinson and Smith, 2001). Diffusion tensors at each voxel were fitted by the algorithm included in the FMRIB's Diffusion Toolbox of FMRIB's Software Library [FSL v. 5.0 (Smith et al., 2004)]. Non-brain parts were removed with the brain extraction tool (BET) (Smith, 2002). FA, MD, and AD and RD to the principal diffusion direction were computed for the whole brain.

Tract-based spatial statistics method (Smith et al., 2007) was used to search for the white matter alterations. All subjects' FA images were aligned into a common space, using the non-linear registration tool, FNIRT, which uses a b-spline representation of the registration warp field. A mean FA image was created and the threshold set at FA = 0.2, deriving a mean FA skeleton that represented the centers of all tracts common to the group. Each subjects' aligned FA data were then projected onto this skeleton and the resulting data fed into voxel-wise cross-subject statistics. Modeling and inference using standard GLM design set-up was accomplished using permutation-based cluster analysis (n = 5000) as implemented in the FSL software package (Nichols and Holmes, 2002). The regressors of the GLM analysis coded for group membership or clinical variables. Three consecutive analyses were conducted for group comparisons: MwA vs. controls, MwoA vs. controls and MwA vs. MwoA. Correlation analysis was conducted between diffusivity parameters and disease duration and attack frequency. The regressors (disease duration and attack frequency) were demeaned. With the GLM design negative and positive correlation were calculated. Statistical thresholding was carried out with Threshold Free Cluster Enhancing (Smith and Nichols, 2009).

p-Values less than 0.05 corrected was accepted as significant result. Statistical maps were thresholded for 0.05 and the Johns Hopkins University white-matter atlas was used to identify the anatomical locations of altered regions.

#### TABLE 1 | Demographic and clinical data.

fnana-11-00138 January 13, 2018 Time: 15:58 # 3


#### RESULTS

#### Clinical Variables

The clinical and demographic variables of the patients and the control group are summarized in **Table 1**. No significant differences were found between the age or gender distribution of the groups (p > 0.05). The patients' groups didn't differ in disease duration (p > 0.05) or attack frequency (p > 0.05).

#### Group Differences in White Matter Microstructure MwA vs. Controls

Decreased RD (p < 0.036; peak Z score = 4.545) was found in MwA compared to controls in the corpus callosum, bilaterally in the parieto-occipital and in the cingular white matter. There was a trend showing decreased MD (p < 0.068; peak Z score = 3.541) in MwA in overlapping areas and increased FA (p < 0.061) in the corpus callosum (**Figure 1**).

#### MwoA vs. Controls

No significant alterations were found between MwoA and the control group.

#### MwA vs. MwoA

Higher FA (p < 0.048; peak Z score = 3.974) was found in MwA in the left parieto-occipital white matter (**Figure 2**).

FIGURE 2 | White matter alterations in MwA compared to MwoA. Tract-TBSS indicate increased FA in MwA compared to migraine without. The color bar shows the z-scores of the corrected p-values. The boxplot shows the MD parameters depicted from the affected areas, the central mark is the median and the boxes represent the 25 and 75% percentiles.

### Correlation of Clinical Variables with Diffusion Parameters

In MwA, the AD negatively correlated with disease duration in the left superior longitudinal fascicle (left parieto-occipital region) and in the left corticospinal tract (p < 0.036; peak Z score = 5.765). The estimated lifetime attack number showed a negative correlation to the AD in the right superior longitudinal fascicle (p < 0.048; peak Z score = 5.621) in the MwA group (**Figure 3**). There was no significant correlation of the clinical variables in the MwoA group.

### DISCUSSION

In this study we provided evidence for interictal white matter microstructural alterations in MwA, a difference not appearing in MwoA patients.

In our very recent study, the amplitude of the resting state functional MRI activity fluctuation was higher in MwA compared to MwoA (Farago et al., 2017). One possible source of the difference between the two subgroups of migraine patients could be the presence of CSD in MwA. Functional imaging studies showed the signature of the slow depolarization wave spreading over the cortex during aura phase (Hadjikhani et al., 2001). No

fnana-11-00138 January 13, 2018 Time: 15:58 # 4

similar sign of slow depolarization was detected in MwoA patients. CSD induces the release of neurotransmitters that leads to neuroinflammation, glial cell activation (Charles and Brennan, 2009) and oxidative stress (Shatillo et al., 2013). Moreover, there is evidence that CSD modifies circulation in the brain, altering its susceptibility to ischaemia (Ayata and Lauritzen, 2015). All these processes might result in white matter abnormalities. However, while molecular markers of cell damage such as S100B and neuron-specific enolase can be detected in the serum of migraine patients (Teepker et al., 2009; Yilmaz et al., 2011), no difference was found between MwA and MwoA patients that could reflect the structural damage that we have found.

Based on the typical clinical course of the symptoms during the aura, it is thought that the CSD originates from the occipital cortex and spreads over the most of the postcentral brain (Petrusic and Zidverc-Trajkovic, 2014). In our study, the signature of altered microstructure was found in the inferior fronto-occipital fasciculus, thalamo-cortical tracts and the corpus callosum of MwA as compared to controls. Differences between MwA and MwoA patients were detected in the occipito-parietal region that might be coherent with the proposed spread of CSD during the aura phase. Although, it cannot fully explain the topical nature of the differences seen between MwA and MwoA. CSD in MwoA may affect silent areas in the brain and therefore no aura symptoms develop (Ayata, 2010; Cosentino et al., 2014). Furthermore, we cannot exclude that the detected structural alterations are related to the fact that CSD is causing robust changes to brain circulation (Lauritzen, 1987).

Another possible option behind the microstructural alterations might be the cortical hyperexcitability in migraine (Gawel et al., 1983; Aurora et al., 1998; Afra et al., 2000; Antal et al., 2005; Chadaide et al., 2007; Pierelli et al., 2013). As was shown in a recent meta-analysis, the TMS-evoked phosphene threshold is lower in MwA and the prevalence of phosphenes was also higher in MwA (Brigo et al., 2012). Cross-modal perception in MwoA and MwA suggests abnormal cortical excitability, fnana-11-00138 January 13, 2018 Time: 15:58 # 5

which is more affected in MwA patients (Brighina et al., 2015). Recent reports also pointed out that the visual evoked potential measured hyperexcitability predominantly true for MwA (Sand et al., 2008; Coppola et al., 2015).

It is disputed whether there is an alteration in white matter microstructure in MwA and MwoA. In a recent study, Tedeschi et al. (2016) detected no difference in diffusivity parameters in MwA and MwoA. An earlier study by Granziera et al. (2006) found no difference between MwA and MwoA groups using voxel-based morphometry-style analysis of diffusion parameters. In our previous study in a mixed group of episodic MwA and MwoA, we presented a relatively smaller pre-frontal white matter alteration (Szabo et al., 2012). Similarly to our study, using TBSS approach, diffusivity parameter changes were published in MwoA (Yuan et al., 2012; Yu et al., 2013). Neeb et al. (2015) found no DTI changes in chronic and episodic MwoA. These studies are consistent in finding reduced FA in the white matter of patients with MwoA. In our cohort, there was no significant difference in the white matter diffusion parameters of MwoA and controls, like in Tedeschi's paper (Tedeschi et al., 2016). However, the MwA patients showed a significantly higher FA in extensive white matter regions, which differs from that study. The importance of the number of the diffusivity directions and how the tensors are calculated is underlined by the literature (Barrio-Arranz et al., 2015). We used the same approach to analyze the DTI data, but Tedeschi's group measured DTI with 32-diffusivity direction once and in our study, we used 60 directions DTI measured twice which might account for the different results. Additionally, the variant MRI parameters (Landman et al., 2007) used in clinical studies and the number of gradient directions have measurable effects on estimated parameters (Barrio-Arranz et al., 2015).

In our previous paper we proposed two alternative hypotheses to explain the white matter microstructural changes: (i) repeated painful conditions or increased cortical excitability might cause maladaptive plastic changes or alternatively (ii) cortical hyperexcitability and CSD might cause degenerative changes through glutamatergic excitotoxicity and mitochondrial dysfunction (Longoni and Ferrarese, 2006; Moskowitz, 2007; Sas et al., 2010). It seems that in MwA patients the first alternative is dominant, represented by increased FA. Use-dependent plasticity-related gray matter morphological alterations, together with the white matter microstructural changes, were reported earlier (Draganski et al., 2004; Boyke et al., 2008; Teutsch et al., 2008; Scholz et al., 2009; Lerch et al., 2011). The underlying mechanism is still disputed. In learning, where repeated stimuli as "mind training" occur, it is thought to be related to the locally enhanced myelination that is represented by an increase in FA (Wedeen et al., 2005; Blumenfeld-Katzir et al., 2011; Lerch et al., 2011; Sampaio-Baptista et al., 2013). Recently, it was suggested that a correlation exists between nerve conduction velocity and locally measured FA (Akhtari et al., 2006; Brienza et al., 2014; Heckel et al., 2015; Wang et al., 2015). We hypothesize that repeated painful attacks or cortical hyperexcitability induce used-dependent plasticity in MwA patients, which in the white matter is represented by increased FA, a signature of a more compact structure.

The longer the history of migraine, the lower the white matter AD. This finding is similar to that described by Yu et al. (2013). Altered AD was reported in axonal loss (Pierpaoli et al., 2001; Kim et al., 2006; Sun et al., 2006) and molecular markers also point to neuronal and glial damage in migraine (Yilmaz et al., 2011). A negative correlation between disease duration and AD might be the sign of chronic myelin and axonal damage (Kincses et al., 2013; Benitez et al., 2014; Gregory et al., 2015), such as in neurodegenerative disorders. Another explanation might be that migraine pathology doesn't affect the brain areas equally and some pathways become more robust due to plasticity.

### CONCLUSION

Migraine is a heterogeneous disease and our results indicate that degenerative and maladaptive plasticity coexist in the disease. The variable diffusion profiles described in the current study and earlier investigations are a consequence of these processes. Our results also point to a possible difference in the pathomechanism of MwA and MwoA, suggesting a separate investigation of these two forms.

### AUTHOR CONTRIBUTIONS

LV, ZTK, and NS planned the project and formulated the study hypothesis. ÁP, JT, and DS recruited the patients. BK, GC, ET, and KK organized and carried out the MRI measurements. BT, DV, and NS collected the clinical data. AK, PF, and NS analyzed the MRI data and carried out the statistical analysis. NS, DV, ZTK, ÁP, DS, and JT formalized the discussion of the results. NS, ZTK, DV, LV, and ÁP wrote the manuscript.

### ACKNOWLEDGMENTS

The study was supported by the "MTA-SZTE Neuroscience Research Group" (Grant No. GINOP-2.3.2-15-2016-00034), EU-funded Hungarian (Grant No. EFOP-3.6.1-16-2016-00008) and the National Brain Research Program 2 (Grant No. 2017- 1.2.1-NKP-2017-00002). NS was supported by the Bolyai Scholarship Program of the Hungarian Academy of Sciences. We acknowledge the core facility MAFIL of CEITEC supported by the MEYS CR (LM2015062 Czech-BioImaging). AK, PF, and ET were supported by the UNKP-17-3 New National Excellence Program of the Ministry of Human Capacities. This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 734718 (CoBeN).

## REFERENCES

fnana-11-00138 January 13, 2018 Time: 15:58 # 6


morphometry studies. Neuroscience 299, 88–96. doi: 10.1016/j.neuroscience. 2015.04.066


fnana-11-00138 January 13, 2018 Time: 15:58 # 7


**Conflict of Interest Statement:** 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.

Copyright © 2018 Szabó, Faragó, Király, Veréb, Csete, Tóth, Kocsis, Kincses, Tuka, Párdutz, Szok, Tajti, Vécsei and Kincses. 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) or licensor 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.

# Human Thalamic-Prefrontal Peduncle Connectivity Revealed by Diffusion Spectrum Imaging Fiber Tracking

Chuanqi Sun<sup>1</sup> , Yibao Wang<sup>2</sup> , Run Cui <sup>2</sup> , Chong Wu<sup>3</sup> , Xinguo Li <sup>2</sup> , Yue Bao<sup>2</sup> and Yong Wang<sup>2</sup> \*

<sup>1</sup> Department of Innovative Medical Photonics, Preeminent Medical Photonics Education & Research Center, Institute for Photonics Research, Hamamatsu University School of Medicine, Hamamatsu, Japan, <sup>2</sup> Department of Neurosurgery, First Affiliated Hospital of China Medical University, Shenyang, China, <sup>3</sup> Department of Neurosurgery, Anshan Central Hospital, Anshan, China

The thalamic-prefrontal peduncle (TPP) is a large bundle connecting the thalamus and prefrontal cortex. The definitive structure and function of the TPP are still controversial. To investigate the connectivity and segmentation patterns of the TPP, we employed diffusion spectrum imaging with generalized q-sampling reconstruction to perform both subject-specific and template-based analyses. Our results confirmed the trajectory and spatial relationship of the TPP in the human brain and identified the connection areas in the prefrontal cortex. The TPP-connecting areas identified based on Brodmann areas (BAs) were BAs 8–11 and 45–47. Based on the automated anatomical atlas, these areas were the medial superior frontal gyrus, superior frontal gyrus, middle frontal gyrus, pars triangularis, pars orbitalis, anterior orbital gyrus, and lateral orbital gyrus. In addition, we identified the TPP connection areas in the thalamus, including the anterior and medial nuclei, and the lateral dorsal/lateral posterior nuclei. TPP fibers connected the thalamus with the ipsilateral prefrontal BAs 11, 47, 10, 46, 45, 9, and 8 seriatim from medial to lateral, layer by layer. Our results provide further details of the thalamic-prefrontal peduncle structure, and may aid future studies and a better understanding of the functional roles of the TPP in the human brain.

#### Edited by:

Laurent Petit, Centre National de la Recherche Scientifique (CNRS), France

#### Reviewed by:

Ilyess Zemmoura, François Rabelais University, France Juan C. Fernandez-Miranda, University of Pittsburgh Medical Center, United States

#### \*Correspondence:

Yong Wang wangyongdl@126.com

Received: 19 September 2017 Accepted: 15 March 2018 Published: 17 April 2018

#### Citation:

Sun C, Wang Y, Cui R, Wu C, Li X, Bao Y and Wang Y (2018) Human Thalamic-Prefrontal Peduncle Connectivity Revealed by Diffusion Spectrum Imaging Fiber Tracking. Front. Neuroanat. 12:24. doi: 10.3389/fnana.2018.00024 Keywords: thalamic-prefrontal peduncle, diffusion spectrum imaging, fiber tractography, connectivity pattern, white matter

### INTRODUCTION

Thalamic peduncles, or radiations, are a large bundle of fan-like fibers connecting the thalamus and cortex, and well-described in both human and monkey brains. Based on the connected cortical area, thalamic peduncles are commonly divided into five components: the anterior, superior, lateral (or extracapsular), posterior, and inferior peduncles (de Castro et al., 2005; Schmahmann and Pandya, 2006; Serra et al., 2017). A classic neuroanatomy study using white matter fiber microdissection defined the fibers originating in the thalamus and radiating along the anterior limb of the internal capsule toward the prefrontal cortex as the anterior thalamic peduncle (ATP), and the fibers radiating from the thalamic pulvinar to the tip of the temporal lobe as the inferior thalamic peduncle (ITP) (Serra et al., 2017). However, a recent autoradiography study in monkeys, as well as human studies, reported that the ITP also contains a large bundle of fibers extending to the orbital cortex (Schmahmann and Pandya, 2006; Jiménez et al., 2007, 2013). To avoid confusion, we

consider all fibers connecting the thalamus and prefrontal cortex as one fiber bundle, termed the thalamic-prefrontal peduncle (TPP) and composed of the ATP and parts of the ITP.

In recent years, diffusion imaging techniques, such as diffusion tensor imaging (DTI) and improved high-definition fiber tracking (HDFT), have provided non-invasive methods for investigating brain white matter fibers and connectivity in vivo (Basser et al., 2000; Fernandez-Miranda et al., 2012). Previous DTI studies have delineated the trajectories of the thalamic peduncles. The fibers originate from the thalamus and project to the brain cortex through the internal capsule. Although DTI reveals fiber trajectory, it does not provide detailed information on the connectivity patterns and precise termination areas of the thalamic peduncles (Wakana et al., 2004). An HDFT study tracking white matter fibers and resolving complex fiber crossings described the thalamocortical projection system (Fernandez-Miranda et al., 2012). The study observed the thalamic fibers radiating to the frontal, central, parietal, and occipital cortical areas, and identified the origin/termination areas of the fibers in the thalamus: the medial and anterior portions of the thalamus connected to the orbitofrontal cortex, the ventral anterior and lateral potions to the prefrontal cortex, and most of the posterior potion to the parieto-occipital cortex. However, this study still did not reveal the precise connected cortical regions of the fibers (Fernandez-Miranda et al., 2012). Another study, based on probabilistic tractography algorithms, showed similar results to those obtained with HDFT; however, it also failed to identify the detailed cortical connections of the thalamic peduncles. The study treated the entire prefrontal cortex as a single region, without subdivision, and only characterized the connections between the thalamic mediodorsal and ventral anterior nuclei, and the prefrontal cortex (Behrens et al., 2003).

The existence of the TPP has been established in classic neuroanatomy, autoradiography, and diffusion imaging studies of monkey and human brains. However, the detailed TPP structure and connectivity patterns, such as specific terminal areas in the prefrontal cortex, remain unclear. In the present study, we aimed to identify the connectivity patterns of the human TPP. To overcome the limitations of DTI studies (Alexander and Barker, 2005; Le Bihan et al., 2006), we employed a high-resolution form of white matter tractography, diffusion spectrum imaging (DSI) (Wedeen et al., 2005, 2008), reconstructed by generalized q-sampling imaging (GQI) (Yeh et al., 2010; Zhang et al., 2013), to sufficiently resolve the complex fiber crossings involved in the TPP.

### MATERIALS AND METHODS

#### Participants

Nine neurologically healthy adults (4 males; all right-handed; age range 22–34) took part in this experiment. All participants were prescreened prior to testing to rule out any contraindications to MR imaging. Written informed consent was obtained from all participants prior to testing. The procedures used were approved by the institutional review board, including the ethics committee of China Medical University.

#### Imaging Acquisition and Reconstruction

DSI data were acquired on a 3T Tim Trio System (Siemens) with a 32-channel head coil. A head stabilizer was used to prevent head motion. Image acquisition involved a 43 min, 257-direction scan using a twice-refocused spin-echo echo-planar imaging sequence and multiple q-values (echo time [TE] = 157 ms, repetition time [TR] = 9,916 ms, field of view [FoV] = 231 × 231 mm, voxel size = 2.4 × 2.4 × 2.4 mm, bmax = 7000 s/mm<sup>2</sup> ) (Wedeen et al., 2005). We also included high-resolution anatomical imaging for anatomical comparisons, employing a 9 min T1-weighted axial magnetization prepared rapid gradient echo (MPRAGE) sequence (TE = 2.63 ms, TR = 2,110 ms, flip angle = 8 ◦ , 176 slices, voxel size = 0.5 × 0.5 × 1.0 mm, FoV = 256 × 256 mm). DSI data were reconstructed with a generalized Q-sampling imaging approach (Yeh et al., 2010). The orientation distribution functions (ODFs) were reconstructed to 362 discrete sampling directions and a mean diffusion distance of 1.2 mm.

### Human Connectome Projection (HCP)-842 Template

In addition to subject-specific analyses, we conducted fiber tracking on a publicly available high-definition fiber tracking template (HCP-842 atlas). The data were averaged from a total of 842 subjects included in the HCP dataset maintained by the WU-Minn HCP Consortium (2015, 900-subject release) and distributed under the WU-Minn HCP open access data use terms (Van Essen et al., 2013). Diffusion images were acquired using a multishell diffusion scheme. The b-values were 1000, 2000, and 3000 s/mm<sup>2</sup> , with 90 diffusion sampling directions in each case. The in-plane resolution was 1.25 mm, and the slice thickness was 1.25 mm. The diffusion data were reconstructed in the MNI space using generalized q-space reconstruction (Yeh and Tseng, 2011) to obtain the spin distribution function. A diffusion sampling length ratio of 1.25 was used, and the output resolution was 2 mm. The atlas was constructed by averaging the SDFs of 842 individual subjects. The HCP-842 template is available freely for download at http://dsi-studio.labsolver.org.

#### Fiber Tracking and Analysis

The DSI Studio software, an open-source diffusion algorithm analysis tool freely downloadable at http://dsi-studio.labsolver. org, was used for fiber tracking. Rather than adopt a wholebrain fiber-tracking procedure, we chose to use an approach based on ODF-streamlined regions of interest (ROIs) (Yeh et al., 2010). In voxels with multiple fiber orientations, fiber tracking was initiated separately for each orientation, and fiber progression continued with a step size of 1.2 mm (subjects) or 1 mm (template; half the spatial resolution), a minimum fiber length of 20 mm, and a turning angle threshold of 60◦ . If multiple fiber orientations existed in the current progression location, the fiber orientation nearest to the incoming direction and forming a turning angle smaller than 60◦ was selected to determine the next moving direction. To smoothen each track, the next moving directional estimate of each voxel was weighted by 20% of the previous incoming direction and 80% of the nearest fiber orientation (Wang et al., 2013; Fernández-Miranda et al., 2015).This progression was repeated until the quantitative anisotropy (QA) (Yeh et al., 2010) of the fiber orientation dropped below a preset threshold (0.03–0.06 depending on the subject) or until no fiber was selected within the 60◦ angular range in the progression. The tracking was stopped when it had generated 40,000 fiber tracts.

To reconstruct the right and left TPP, we used bilateral prefrontal cortical ROIs, including Brodmann areas (BAs) 8–11 and 44–47, as seeding regions (**Figure 1**). Another ROI was set at the left or right thalamus as an end region. To quantify and compare the right and left sides, as well as the contribution of each seeding region, we counted the numbers of voxels occupied by the fiber trajectories (streamlines). In addition, to understand the spatial relationship of the TPP with adjacent fiber tracts, the inferior frontal-occipital fasciculus (IFOF), anterior part of the corpus callosum (A-CC), and corticospinal tracts (Cort-Sp) were reconstructed using the same fiber tracking parameters as in the TPP reconstruction.

To determine the cortical terminations and segmentation of the TPP, we imported cortical prefrontal areas from two different atlases, the BA and improved automated anatomical (AAL2)

FIGURE 1 | Fiber tractography of the thalamic-prefrontal peduncle (TPP; subject No. 8). (A) Seeding regions of interest (ROIs) were placed on the prefrontal cortex, including Brodmann areas (BAs) 8–11 and 44–47. The thalamus (Tha) was set as the ending ROI. The ROIs were overlaid on the white-matter surface. (B) Axial slice at the thalamus level on a quantitative anisotropy (QA) color map; Ending ROIs were placed at left thalamus (blue) and right thalamus (red). (C) Posterior-superior view of the TPP and ROIs. (D–F) Lateral, superior, and anterior views of the TPP, overlaid on T1 MRI slices.

atlases (Tzourio-Mazoyer et al., 2002; Rolls et al., 2015). DSI Studio used linear transformations to register both atlases on a subject's diffusion space.

#### Statistical Analysis

Statistical analyses were carried out using SPSS 22.0 (SPSS, Chicago, IL). T-tests were used to evaluate parametric differences. P ≤ 0.05 was considered to be statistically significant. Continuous variables are presented as means ± standard deviation.

### RESULTS

### Human TPP Trajectory

Subject-specific fiber tractography consistently showed a large bundle of fibers projecting anteriorly from the thalamus toward the prefrontal lobe, which included most of the selected seeding regions (**Figures 1C–F**). Template (HCP-842) tractography showed similar results (**Figure 2**). This distinct group of fibers originated from the thalamus and projected along the anterior limb of the internal capsule to the prefrontal cortex, primarily including the superior frontal gyrus (SFG), orbital-frontal cortex (OFC), and inferior frontal gyrus (IFG). A qualitative study of 20 hemispheres (18 hemispheres from 9 subjects and 2 hemispheres from the HCP-842 template) showed nearly identical trajectories (**Figure 3**).

### Spatial Relationship of the TPP With Adjacent Fiber Tracts and Structures

To visualize the anatomical spatial relationship of the TPP with adjacent tracts, the Cort-Sp, IFOF, and A-CC were reconstructed using in vivo fiber tractography.

The Cort-Sp has been well-described and is situated superficially to the TPP. In our analyses, the main trunk of the TPP was located medially to the Cort-Sp tract and projected anteriorly to the prefrontal cortex (**Figures 4A–C**).The IFOF is a large bundle of fibers that tracks from the prefrontal region to the occipital region via the ventral part of the external/extreme capsule. The prefrontal part of the IFOF was situated immediately ventrally and superficially to the prefrontal part of the TPP fibers (**Figure 4A**). Importantly, the interconnected prefrontal cortex regions of the TPP and IFOF overlapped partially at the IFG and orbital gyrus. The corpus callosum (CC) is a large bundle of commissural fibers connecting the left and right brain hemispheres. The A-CC was located medially to the TPP (**Figures 4B,C**).

The Free Surfer parcellations (Desikan et al., 2006) of the caudate, putamen, and pallidum nuclei were overlaid on the TPP

FIGURE 3 | Comparison of the thalamic-prefrontal peduncles (TPPs) between 9 subjects (subject: 1–9) and the HCP-842 template. Both the left (Left lateral view) and right (Right lateral view) TPP of each subject are presented. Similar trajectories were observed in 20 hemispheres (18 hemispheres from 9 subjects; 2 hemispheres from the HCP-842 template).

FIGURE 4 | Spatial relationship of the thalamic-prefrontal peduncle (TPP) with adjacent fiber tracts and structures. (A) Lateral view shows the spatial relationships of the left TPP (red) with the inferior frontal-occipital fasciculus (IFOF; green) and corticospinal tracts (Cort-Sp; blue). (B) Posterior-medial view of the left TPP, Cort-Sp, and anterior corpus callosum (A-CC; yellow). (C) Inferior view of the TPP, Cort-Sp, and A-CC. (D) Inferior view of the bilateral TPP shows the fiber bundles projecting to the prefrontal cortex through the anterior limb of the internal capsule. (E) Inferior view reveals the spatial relationships of the TPP with the thalamus (Tha), and caudate (Cau), putamen (Put), and pallidum (Pal) nuclei. (F) Posterior view of the bilateral TPP.

TPP-connected areas in the prefrontal cortex included Brodmann areas (BAs) 8, 9, 10, 11, 45, 46, and 47. (B) Segmentation of the TPP fibers based on prefrontal BAs. Different fibers are assigned different colors. (C–I) Views of each tract based on the segmentation. Fibers of each segmentation are colored according to (A,B). (**Figure 4E**). Our results showed the TPP pathway to be adjacent to the head of the caudate nucleus, the putamen, and the pallidum nucleus. The fiber bundles projected to the prefrontal cortex through the anterior limb of the internal capsule (**Figures 4C–F**).

#### TPP Segmentation and Connectivity

To investigate the connections of the TPP fibers between the prefrontal cortex and thalamus, we identified different prefrontal cortical regions based on BAs and the AAL2 atlas. The areas connected to the prefrontal cortex primarily included BAs 8, 9, 10, 11, 45, 46, and 47 (**Figure 5**). TPP fibers terminated bilaterally at BAs 9, 10, and 11 in all 9 subjects. Fibers converged bilaterally at BAs 8, 45, and 46 in 8 out of 9 subjects. Fibers terminated at BA 47 on the left side in 7, and on the right side in 6 out of 9 subjects. Additionally, some fibers terminated at BA 44 on the left side in only 2, and on the right side in only 1 out of 9 subjects. Imaging using the HCP-842 template showed similar results (**Table 1**, **Figure 5**). TPP segmentation based on BAs is shown in **Figure 5**.

We used the AAL2 atlas to observe TPP segmentation and connectivity patterns as an alternative approach to understanding the anatomical regions of TPP fiber termination (**Figure 6**). Seven different TPP termination areas in the prefrontal cortical region were identified: the SFG, medial SFG (MSFG), middle frontal gyrus (MFG), pars triangularis (PTri), pars orbitalis (POrb), anterior orbital gyrus (AOG), and lateral orbital gyrus (LOG). The anterior view of the TPP endpoints has the appearance of the letter "L" or a "fish hook" (**Figures 2C,D**, **6B**). A small bundle of fibers terminated at the caudal MSFG, and a large bundle connected with the middle and rostral SFG. The TPP endpoints were located laterally to the OFC (i.e., AOG and LOG), and in the IFG (i.e., POrb and PTri) and rostral MFG (**Figure 6**).

To characterize the TPP origin/termination areas in the thalamus, the TPP fibers and their endpoints were overlaid with


The tract (and sub-tract) volumes were counted by the number of voxels occupied by the fiber trajectories (streamlines) and are reported in milliliters (ml). The figures in parentheses represent the relative volume of each tract and subtract in percentages (%), where the left and right TPP together represent 100% of the volume. BA, Brodmann Area. "—" indicates that the mean of the tract volumes of BA 44 were not calculated because only three out of 18 hemispheres showed connected tracts.

FIGURE 6 | Cortical connectivity and segmentation of the thalamic-prefrontal peduncle (TPP) based on the AAL2 atlas (HCP-842 template). (A) The TPP-connected areas in the prefrontal cortex included the SFG, MSFG, MFG, Ptri, Porb, LOG, and AOG. (B) Anterior-lateral view: Endpoints of the TPP were overlaid on the white-matter surface. (C–I) Views of each tract based on the segmentation. Fibers of each segmentation are colored according to (A). MSFG, medial superior frontal gyrus; SFG, superior frontal gyrus; MFG, middle frontal gyrus; PTri, pars triangularis; POrb, pars orbitalis; AOG, anterior orbital gyrus; LOG, and lateral orbital gyrus.

thalamic T1 MRI slices (**Figure 7**). The origin/termination areas in the thalamus were separated into three parts: anterior, ventral, and dorsal (**Figure 7B**). These parts correspond to the medial (MD) and anterior (Ant) nuclei, ventral-posterior part of the MD, and lateral dorsal/lateral posterior nuclei (LD/LP), respectively (**Figures 7D–F**). To analyze the connectivity patterns of the human TPP fibers, fibers were separated into several subtracts based on the prefronta cortical BAs involved (**Figure 8**). The left TPP endpoints showed TPP tract connections with the prefrontal cortex and thalamus layer by layer. The location sequence of the subtracts in the thalamus was as follows (from medial to lateral): BA 11, BA 47, BA10, BA 46, BA 45, BA9, and BA8 (**Figures 8B,C**). The BA 8 tract was located at the most lateral part of the TPP, the BA 9 tract ran medially to the BA 8 tract, and the BA 45 tract was located medial-ventrally to the BA9 tract. The BA 46 tract ran medially to the BA 45 tract and split into two parts (dorsal and ventral). The BA 10 tract ran medially to the BA 46 tract. The BA 47 tract ran between the BA 11 and BA 10 tracts, turning laterally to the cortex-crossing part of the BA 10 tract (**Figures 8D–F**).

## Quantitative Study of the TPP

The TPP volumes in 9 subjects and the HCP-842 template are presented in **Table 1**. The results were similar between the subject-specific and template-based approaches. Statistical analysis of the subject-specific approach revealed a significant difference between the left and right TPP in tract volume (P = 0.03). However, when subcomponents of the TPP were analyzed, only BA 9 showed a significant difference between the left and right hemispheres (P = 0.048). The other subcomponents showed no significant differences between the left and right hemispheres (P > 0.05) (**Table 2**).

#### DISCUSSION

Previous studies have delineated the TPP fibers in monkey and human brains using different methods, including fiber microdissection, autoradiography, and diffusion imaging studies. However, these reports do not contain detailed information on the connections between the thalamus and prefrontal

thalamus corresponds to the medial (MD) and anterior (Ant) nuclei of the thalamus. (E) The ventral origin/termination area corresponds to the ventral-posterior part of

the MD. (F) The dorsal origin/termination area of TPP fibers in the thalamus corresponds to the lateral dorsal/lateral posterior nuclei (LD/LP).

cortex. Here, we further characterize the human TPP in vivo using the DSI technique. We also report the fiber trajectories, spatial relationships, cortical connectivity, segmentation, and a quantitative assessment of the human TPP.

We analyzed the TPP using two approaches, subject-specific and template-based, to enhance the reliability of our results. The TPP connections presented here were first observed in 9 subjects (18 hemispheres) and confirmed by an HCP-842 template, the average of 842 subjects' HCP data. Using this template, we obtained results from a large sample simultaneously. Unlike the subject-specific approach, the template approach does not allow for the assessment of individual differences. However, averaging data from a larger sample does provide unique information on TPP connectivity patterns. Significantly, the two approaches revealed similar connectivity patterns in the present study.

Descriptions of the thalamic peduncles date back to a century ago, when Charles Judson Herrick first delineated the fibers connecting the thalamus and cortex (Herrick, 1915). Autoradiographic technology, the "gold standard" of imaging, enabled the study of white matter fiber connectivity. An autoradiography study of the monkey brain described in detail the thalamic peduncles, and showed that all regions of the prefrontal cortex sent fibers to the thalamus via the anterior limb of the internal capsule to form the ATP, whereas fibers from the orbital cortex connected to the thalamus through the ITP (Schmahmann and Pandya, 2006). However, recent evidence, such as a report on the middle longitudinal fascicle (Wang et al., 2013), appeared to show that the white matter fiber pathways in humans and monkeys were partially different. Because the present study was based on the human brain in vivo, our results are likely to be more reliable. Previous fiber microdissection studies have identified the TPP trajectory directly (de Castro et al., 2005; Serra et al., 2017), showing that the TPP fibers originating from the thalamus radiate anteriorly and constitute part of the anterior limb of the internal capsule. Studies comparing DSI with classic fiber microdissection and autoradiography techniques have shown DSI to be in good agreement with the other methods both in human and monkey brain white matter (Schmahmann et al., 2007; Fernandez-Miranda et al., 2012; Wang et al., 2013; Fernández-Miranda et al., 2015). However, because the thalamocortical fibers, which join into the internal capsule, are difficult to distinguish from the other fibers by microdissection, the fiber dissection studies may not precisely describe the continuous fiber trajectories from the thalamus to brain cortex. Because we utilized the DSI technique, which can sufficiently resolve the complex fiber crossings, we were able to track the fibers from the thalamus to prefrontal cortex continuously, and to identify the termination areas of the TPP fibers in the prefrontal cortex. Our results show that the TPP fibers project anteriorly from the thalamus toward the prefrontal lobe along the anterior limb of the internal capsule. The connected areas in the prefrontal cortex are BAs 8, 9, 10, 11,

only described fan-like fibers projecting to the SFG. Because of the influence of the crossing fibers, such as in the A-CC, TPP fiber trajectories tracked by DTI may not include all cortical connections. Because DSI represents a better method of fiber tracking, we were able to characterize the spatial relationship of the TPP with adjacent fiber tracts (Cort-Sp, IFOF, and A-CC) (**Figure 4**). HDFT studies have shown the origins/terminations of the thalamic peduncles: the fibers originating from the orbitalfrontal region are linked to the most medial and anterior portions

Pairwise difference t df P

Lower limit Upper limit

TABLE 2 | The statistical analysis data of the volume of the TPP and its subcomponents.


medial-ventrally to the BA9 tract. (E) The BA 46 tract runs medially to the BA 45 tract and splits into two parts (dorsal and ventral tracts). (F) The BA 10 tract runs medially to the BA 46 tract. (G) The BA 47 tract runs between the BA 11 and BA 10 tracts, then turns laterally to the cortex-crossing part of the BA 10 tract.

L—R 5.26 5.98 1.99 0.66 9.86 2.64 8 0.030 L-BA 8—R-BA 8 −0.42 1.28 0.43 −1.41 0.56 −0.99 8 0.349 L-BA 9—R-BA 9 0.89 1.14 0.38 0.01 1.76 2.34 8 0.048 L-BA 10—R-BA 10 0.63 3.02 1.01 −1.69 2.95 0.62 8 0.551 L-BA 11—R-BA 11 0.63 3.30 1.10 −1.91 3.17 0.57 8 0.582 L-BA 45—R-BA 45 1.75 3.11 1.04 −0.64 4.14 1.69 8 0.129 L-BA 46—R-BA 46 0.62 2.64 0.88 −1.41 2.65 0.71 8 0.498 L-BA 47—R-BA 47 1.07 2.88 0.96 −1.15 3.28 1.11 8 0.299

Mean value SD Standard error 95% confidence interval

of the thalamus, those from the prefrontal area to the ventral anterior and lateral regions of the thalamus, central lobule fibers to the ventral posterior zone, and parietal-occipital fibers to the posterior and inferior zones of the thalamus (Fernandez-Miranda et al., 2012). Our study describes the origin/termination areas of the TPP based on the thalamic structure. The TPP links to the Ant, MD, and LD/LP portions of the thalamus, similar to the HDFT studies and in agreement with neurohistological studies (Nieuwenhuys et al., 2008; Fernandez-Miranda et al., 2012). While DTI, HDFT, and DSI utilize deterministic tractography algorithms, a study of the connections between the human thalamus and cortex based on a probabilistic approach showed that the fibers linked the prefrontal cortex with the MD, ventral anterior (VA), and parts of the anterior portion of the thalamus (Behrens et al., 2003), largely in agreement with our results. A functional MRI study also showed functional connectivity between the prefrontal cortex and mediodorsal and anterior nuclear areas of the thalamus, similar to our results (Zhang et al., 2008).

Although the thalamic peduncles have been characterized in detail, their precise cortical connections with the thalamus have remained unclear. For comparison and a better understanding of the connections, we used two prefrontal cortical maps (from the BA and AAL atlases) to describe the precise TPP connections with the region and discovered a novel TPP connectivity pattern. Our results show that the TPP fibers connect the prefrontal cortex with the thalamus layer by layer based on the prefrontal BA atlas. The tract which connects the orbital-frontal area (BA11) with the thalamus is located at the most medial portion of the TPP, the SFG tract (BA8, 9) at the most lateral part of the TPP, and the tracts projecting to the IFG at the middle of the TPP (**Figure 8**). These results were obtained using both the subject-specific and template approaches.

Studies of thalamic function indicate that the thalamus is a relay center for sensory and motor information, awareness, attention, and other cognitive processes such as memory and language (Herrero et al., 2002; de Bourbon-Teles et al., 2014; Dalrymple-Alford et al., 2015). The function of the prefrontal cortex has been researched extensively and includes cognitive abilities, social emotion, executive functioning, motor control, and language (Teffer and Semendeferi, 2012; Nestor et al., 2013). Because the TPP is a large bundle of fibers that connects two important brain regions, the thalamus and prefrontal cortex, it may play a significant role in the functions mentioned above. Interestingly, we found TPP fibers extending to BA45 (PTri), which corresponds to the anterior portion of Broca's area, commonly considered to be the substrate of functional language processing. The brain cortical functions involve white matter

#### REFERENCES


fibers; thus, our findings may provide strong new evidence for the role of the thalamus in language function. In addition to BA45, other large bundles of fibers terminated in the OFC, SFG, and BA47. However, only a small fiber bundle terminated at the very anterior portion of the MFG. Given this pattern of TPP connectivity, we can focus our attention on the correlated areas in future investigations of TPP fiber functions and the relationship between the thalamus and cortex.

Although the DSI technique we employed is better able to assess complex fibers, such as crossing fibers, than is DTI, it still has limitations. For example, when analyzing the TPP, we found some false continuations of fibers at the lateral edge of the TPP; these fibers did not terminate in the thalamus, which should be part of the Cort-Sp. Our study only focused on the connectivity between the prefrontal cortex and thalamus, we may loss some fibers information, especially in the region near the parietal cortex. However, these limitations did not influence the TPP connectivity patterns eventually identified. Previous studies comparing the results of DSI analyses and those of fiber dissections have shown that the former produces more reliable evidence (Fernandez-Miranda et al., 2012; Wang et al., 2013; Fernández-Miranda et al., 2015). Therefore, the results of the present study are likely to be reliable and may provide a foundation for future studies of the thalamic peduncles.

### CONCLUSIONS

We studied the TPP connectivity patterns in the human brain in vivo using both subject-specific and template-based approaches. Detailed connections between the thalamus and prefrontal cortex were described. Our results may further the study of the functions of the thalamus and prefrontal cortex based on fiber connectivity patterns.

### AUTHOR CONTRIBUTIONS

Conceived and designed the experiments: YiW and CS. Performed the experiments: CS, CW, YB, and YoW. Statistical analysis: RC and XL. Data interpretation and figure preparation: CS and YiW. Wrote and revised the manuscript: CS, YoW, and YiW.

### ACKNOWLEDGMENTS

We thank Dr. Seiji Yamamoto for the supporting and helping in the research. This study was supported by grants from the Liaoning Provincial Natural Science Foundation of China (NO. 2014021066 to YoW).


**Conflict of Interest Statement:** 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.

Copyright © 2018 Sun, Wang, Cui, Wu, Li, Bao and Wang. 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 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.

# Awakening Neuropsychiatric Research Into the Stria Medullaris: Development of a Diffusion-Weighted Imaging Tractography Protocol of This Key Limbic Structure

Darren W. Roddy <sup>1</sup> \*, Elena Roman<sup>1</sup> , Shane Rooney <sup>1</sup> , Sinaoife Andrews <sup>1</sup> , Chloe Farrell <sup>1</sup> , Kelly Doolin<sup>1</sup> , Kirk J. Levins <sup>2</sup> , Leonardo Tozzi <sup>3</sup> , Paul Tierney <sup>4</sup> , Denis Barry <sup>4</sup> , Thomas Frodl <sup>3</sup> , Veronica O'Keane<sup>1</sup> and Erik O'Hanlon<sup>1</sup>

<sup>1</sup>REDEEM Group, Department of Psychiatry, Trinity College Dublin, Dublin, Ireland, <sup>2</sup>Department of Anaesthesia, Intensive Care and Pain Medicine, St. Vincent's Hospital, Dublin, Ireland, <sup>3</sup>Department of Psychiatry and Psychotherapy, University of Magdeburg, Magdeburg, Germany, <sup>4</sup>Department of Anatomy, Trinity College Dublin, Dublin, Ireland

#### Edited by:

Silvio Sarubbo, Ospedale Santa Chiara, Italy

#### Reviewed by:

Chiara Maffei, University of Trento, Italy Francesco Sammartino, The Ohio State University, United States

> \*Correspondence: Darren W. Roddy dwroddy@tcd.ie

Received: 31 December 2017 Accepted: 25 April 2018 Published: 08 May 2018

#### Citation:

Roddy DW, Roman E, Rooney S, Andrews S, Farrell C, Doolin K, Levins KJ, Tozzi L, Tierney P, Barry D, Frodl T, O'Keane V and O'Hanlon E (2018) Awakening Neuropsychiatric Research Into the Stria Medullaris: Development of a Diffusion-Weighted Imaging Tractography Protocol of This Key Limbic Structure. Front. Neuroanat. 12:39. doi: 10.3389/fnana.2018.00039 The Stria medullaris (SM) Thalami is a discrete white matter tract that directly connects frontolimbic areas to the habenula, allowing the forebrain to influence midbrain monoaminergic output. Habenular dysfunction has been shown in various neuropsychiatric conditions. However, there exists a paucity of research into the habenula's principal afferent tract, the SM. Diffusion-weighted tractography may provide insights into the properties of the SM in vivo, opening up investigation of this tract in conditions of monoamine dysregulation such as depression, schizophrenia, addiction and pain. We present a reliable method for reconstructing the SM using diffusionweighted imaging, and examine the effects of age and gender on tract diffusion metrics. We also investigate reproducibility of the method through inter-rater comparisons. In consultation with neuroanatomists, a Boolean logic gate protocol was developed for use in ExploreDTI to extract the SM from constrained spherical deconvolution based whole brain tractography. Particular emphasis was placed on the reproducibility of the tract, attention to crossing white matter tract proximity and anatomical consistency of anterior and posterior boundaries. The anterior commissure, pineal gland and mid point of the thalamus were defined as anatomical fixed points used for reconstruction. Fifty subjects were scanned using High Angular Resolution Diffusion Imaging (HARDI; 61 directions, b-value 1500 mm<sup>3</sup> ). Following constrained spherical deconvolution whole brain tractography, two independent raters isolated the SM. Each output was checked, examined and cleaned for extraneous streamlines inconsistent with known anatomy of the tract by the rater and a neuroanatomist. A second neuroanatomist assessed tracts for face validity. The SM was reconstructed with excellent inter-rater reliability for dimensions and diffusion metrics. Gender had no effect on the dimensions or diffusion metrics, however radial diffusivity (RD) showed a positive correlation with age. Reliable identification and quantification of diffusion metrics of the SM invites further exploration of this key habenula linked structure in neuropsychiatric disorders such as depression, anxiety, chronic pain and addiction. The accurate anatomical localization of the SM may also aid preoperative stereotactic localization of the tract for deep brain stimulation (DBS) treatment.

Keywords: stria medullaris, DTI, tractography, CSD, habenula, depression, deep brain stimulation, limbic system

### INTRODUCTION

The dorsal diencephalic conduction system is a pathway that transmits information from the cognitive-emotional forebrain to the regulatory midbrain areas (Sutherland, 1982). This highly conserved system (Beretta et al., 2012) gathers fibers from diverse frontal areas including the septal nuclei (pleasure and motivation), dorsal anterior cingulate (reward-based decision making), lateral hypothalamus (arousal and pain) and basal ganglia (motor and behavioral control). After forming a single unidirectional tract, these fibers project to the habenulae, a set of relay nuclei at the dorso-caudal end of the thalamus (Parent et al., 1981; Geisler and Trimble, 2008). The habenula in turn projects down through the fasciculus retroflexus influencing midbrain monoamines of the raphe nuclei and ventral tegmental area. Through this system the frontolimbic areas can influence midbrain and whole brain monoaminergic tone. The single unidirectional tract connecting the frontobasal areas to the habenula is called the stria medullaris (SM).

The SM (or SM thalami) is a white matter epithalamic tract (**Figure 1**). It emerges just caudal to the anterior commissure as a solitary tract that arches both dorsally and caudally over the inter-thalamic adhesion. It terminates at both lateral and medial habenular nuclei (Parent et al., 1981; Geisler and Trimble, 2008; Ahumada-Galleguillos et al., 2017). The habenular nuclei receive input chiefly through SM (Herkenham and Nauta, 1977). The area of the habenula occupied by fibers that enter the habenula through the SM is considerably larger in humans compared to other species. Post mortem studies suggest that the human SM covers about 30% of the entire cross-sectional area of the habenula, compared to only 12% in rodents (Diaz et al., 2011; Ahumada-Galleguillos et al., 2017). This most likely reflects the increased connectivity of the enlarged human forebrain to the habenula, consistent with its function in relaying information from frontal to midbrain regions.

The habenula outputs either directly or indirectly to midbrain monoaminergic nuclei to affect the release of serotonin (Kalén et al., 1985), dopamine (Gruber et al., 2007) and norepinephrine (Quina et al., 2015). In general, habenula activation reduces monoaminergic output (Boulos et al., 2017). Using human neuroimaging, the habenula has been implicated in aversive decision making (Lawson et al., 2014; Hennigan et al., 2015), reward pathways (Salas et al., 2010), mood (Morris et al., 1999; Ely et al., 2016) and pain (Shelton et al., 2012). Habenula changes have also been found in disorders such as depression (Lawson et al., 2017; Liu et al., 2017) and chronic pain (Erpelding et al., 2014). The unidirectional flow of information from the limbic frontal regions through the SM to the habenula and the apparent direct relationship of habenula activity to neurotransmitter release invite study of the SM in neuropsychiatric disorders such as depression, bipolar disorder, schizophrenia, chronic pain and addiction our (Fakhoury, 2017). Due to its strategic location, the SM may provide information on potential problems upstream to the habenula. However, no study has investigated the role of this key structure in disease. This is due to a lack of techniques to assess this white matter tract in vivo in humans.

Recent advances in diffusion-weighted magnetic resonance imaging (DWI) tractography have allowed identification of white matter tracts with ever-greater precision (Johansen-Berg and Rushworth, 2009; Farquharson et al., 2013). This technique has been particularly useful for identifying large white matter tracts such as the superior longitudinal fasciculus, inferior longitudinal fasciculus and corpus callosum (Mangin et al., 2013). However, parsing out the smaller and more tortuous tracts of the limbic system using tractography can be complicated by methodological difficulties such as resolution issues and crossing fibers (Mori and Aggarwal, 2014). Newer, more advanced technology is aiding with these endeavors. Increasing MR field strengths and the use of faster acquisition sequences can reduce scan times while increasing the signal to noise ratio (Soares et al., 2013). Preprocessing developments such as head motion, distortion and free water correction techniques can also increase the signal to noise ratio (Tournier et al., 2011). Diffusion models such as constrained spherical deconvolution (CSD; Jeurissen et al., 2011), diffusion spectrum imaging (Wedeen et al., 2008) and Q-Ball imaging (Tuch, 2004) are allowing previously hidden crossing, diverging and kissing fibers to be explored. Such developments are facilitating accurate reconstruction of limbic tracts such as the fornix (Christiansen et al., 2017), ventral amygdalofugal pathway (Kamali et al., 2016), stria terminalis (Kamali et al., 2015) and medial forebrain bundle (Hana et al., 2015). Tract diffusion metrics such as fractional anisotropy (FA) and radial diffusivity (RD) may provide information on the condition of limbic white matter tracts in disease (Metzler-Baddeley et al., 2012; Oishi et al., 2012).

The SM as a discrete white matter tract originates just behind the anterior commissure and terminates at the habenula. It does not deviate or branch between these two fixed points. This characteristic makes it a good potential candidate for tractography. However, the SM is particularly challenging to accurately render with DWI tractography. It is short, thin and highly curved (180 degree turn along its course) making it methodologically tricky for DWI tractography to process. It lies deep within the brain making it difficult to isolate from the multiple surrounding white matter structures. Finally, the

as the closest anatomical landmark. As such the protocol routinely cuts the SM slightly short, stopping just before its insertion into the habenula. AC, Anterior

SM and its termination, the habenula, are inconsistently seen on standard T1/T2 MRI due to being ''missed'' between slices (Strotmann et al., 2014; Kochanski et al., 2016).

Commissure; ITA, Inter-thalamic adhesion, SM, Stria Medullaris.

Accurately imaging the SM using tractography could facilitate research into the neuropathology of psychiatric, addiction and pain disorders through the analysis of tract diffusion metrics. Furthermore, its pivotal role in integrating diverse basal forebrain limbic regions into a single tract and relaying directly with the habenula makes the SM a target for deep brain stimulation (DBS) in the treatment of severe depression (Sartorius et al., 2010; Kiening and Sartorius, 2013). Accurately mapping the trajectory may aid the perioperative localization of the tract in stereotactic surgery. This study aimed to reconstruct the SM trajectory in 50 healthy individuals with the limitations of DWI tractography of this region in mind. The signal to noise ratio was maximized by careful preprocessing of the data including signal drift, head motion, eddy current and EPI distortion correction. High Angular Resolution Diffusion Imaging (HARDI) data was acquired and a CSD model was used to create tracts of the entire brain. This algorithm has been shown to overcome some of the known difficulties and shortcomings of the diffusion tensor model and has the ability to model regions with complex fiber orientations such as multiple or crossing fibers within a voxel, allowing more convoluted tracts to be accurately explored. Due to the SM's size and highly curved shape, a high angle threshold and small step size was used. Neuroanatomical expertise was also involved at all stages of protocol design. Reliability between raters was measured to evaluate the consistency of the protocol and age and sex differences of the SM were assessed using diffusion metrics. As such we present here a novel and reliable method for reconstructing the SM in vivo using constrained spherical deconvolution based deterministic tractography.

### MATERIALS AND METHODS

#### Participants

Fifty participants (25 female) were enrolled through an active database of willing participants as part of the REDEEM (Research in Depression: Endocrinology, Epigenetics and neuroiMaging) study at Trinity College Dublin. The mean age and standard deviation (SD) was 29.58 (SD 13.03) with ages of females 29.85 (SD 11.93) and males 29.31 (SD 14.28). Exclusion criteria for study entry included contraindications to MRI, history of illicit substance abuse, head trauma, and any significant medical or psychiatric illness. All participants underwent a Structured Clinical Interview (SCID) for DSM IV, Hamilton Depression Rating Scale (HAM-D; Hamilton, 1960) and Hamilton Anxiety Rating Scale (HAM-A; Hamilton, 1959) administered by a psychiatrist at TCIN to further exclude psychiatric pathology. Subjects were required to have no active or previous SCID diagnosis, a HAM-D score of 8 and below and a HAM-A score of 13 and below.

This study was carried out in accordance with the recommendations of the Tallaght Hospital/St. James Hospital Joint Research Ethics Committee with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by the Tallaght Hospital/St. James Hospital Joint Research Ethics Committee.

#### MRI Acquisition

All data was acquired on a Philips (Best, Netherlands) Intera Achieva 3.0 Tesla MR system (32-channel head coil) at Trinity College Institute of Neuroscience, Dublin. One-hundred and eighty axial high-resolution T1-weighted anatomical images (T1W-IR1150 sequence, TE = 3.8 ms, TR = 8.4 ms, FOV

fornix and extended to the lower thalamus, making sure to encompass the inter-thalamic adhesion if present. The gate should be wide enough to cover most of the "curve" of the thalamus. The inter-thalamic adhesion is labeled. (D) Sagittal and (E) axial-sagittal views of the initial SEED and AND gates. The anatomical fixed points for gate placement: anterior commissure, inter-thalamic adhesion and pineal gland are labeled.

230 mm, 0.898 × 0.898 mm<sup>2</sup> , in-plane resolution, slice thickness 0.9 mm, flip angle alpha = 8◦ ). Whole brain, HARDI (Jones and Leemans, 2011; Jones et al., 2013) was acquired using a spin-echo echo-planar imaging (SE EPI) pulse sequence (TE = 52 ms, TR = 11,260 ms, flip angle alpha = 90◦ ), FOV 224 mm, 60 axial slices, 2 mm<sup>3</sup> isotropic voxels, b-value = 1500 s mm−<sup>2</sup> in 61 non-collinear gradient directions. A single non-diffusionweighted b0 image was also obtained.

#### DWI Preprocessing and Whole Brain Tractography

DWI preprocessing involved the following steps in ExploreDTI<sup>1</sup> (Leemans et al., 2009), a free MATLAB<sup>2</sup> based graphical toolbox for diffusion MRI preprocessing and tractography. Diffusion preprocessing involved the following steps: (a) exporting and standardizing diffusion output files; (b) signal drift correction via linear and quadratic correction (Vos et al., 2017); (c) Gibbs Ringing artifact correction (Perrone et al., 2015); (d)

<sup>1</sup>http://www.exploredti.com

orientation/directionality checks via manual glyph inspection; (e) head motion and eddy current artifact correction using rigid body and affine registration respectively to the non diffusion weighted b0 image (Soares et al., 2013); and (f) EPI deformation correction via affine registration to the T1 image (Irfanoglu et al., 2012). These operations were performed using the latest version toolbox plugins directly within ExploreDTI.

Tractography of the whole brain was generated in ExploreDTI using constrained spherical deconvolution (as a recursive calibration of the response function; Jeurissen et al., 2011). Whole brain seeding was used with seed voxel sizes of 2 mm<sup>3</sup> and a fiber orientation distribution threshold of 0.1. Whole brain tractography was achieved through multiple random seed placements of one seed per voxel (Jones, 2010; Jones and Leemans, 2011). Two separate whole brain tractography files were generated on each subject. First, a sparser one with fewer and less dense streamlines was created to easily scout out the general region of the SM—low fidelity tracts. This was only used for initial gate placement (**Figure 2**) and to identify larger structures that needed to be excluded e.g., the

<sup>2</sup>https://www.mathworks.com/products/matlab.html

fornix. Subsequently, a more complex tract with a higher density of streamlines was created to reconstruct the tract at the latter stages—high fidelity tracts. The step size was set at both 1 mm and 0.5 mm and the maximum angle threshold was set at 45◦ and 89◦ , respectively for low and high fidelity whole brain tracts. All whole brain tract lengths were set between 10–500 mm. The lower threshold was chosen as we assumed the SM would be likely greater than 10 mm (around 30 mm from our own anatomical measurements). The higher threshold was chosen to faithfully reproduce larger known tracts that surround the SM region, such as the fornix (Metzler-Baddeley et al., 2013), stria terminalis (Kamali et al., 2015), thalamocortical radiations (Kubota et al., 2013) and superior thalamic peduncle (Meola et al., 2015). These tracts were subsequently removed with NOT gates. No other tractography stopping criteria apart from angle threshold, FOD threshold and tract length were used.

As mentioned, the low fidelity tracts were used as scout tracts to aid initial gate placement and early cleanup. During our initial investigations into tractography of the SM, we found that having whole brain tracts with a small step size and high angle created overly complex images due to a profuseness of streamlines and multiple interweaving tracts (see **Figure 3** for examples of low and high fidelity tracts prior to cleanup). This created problems when trying to exclude the larger surrounding tracts mentioned above. As a result, we decided that using less complex whole brain tracts with a larger step size and lower angle greatly aided our initial virtual ''dissection'' process. Once the larger known tracts were excluded and we visualized the region where we suspected the SM to be present, we saved

FIGURE 5 | An example of generated tracts post segmentation and splitting (i.e., isolation of the fibers between the anterior and posterior margins only), but prior to neuroanatomical cleanup.

the cleanup gates and applied them to the high fidelity tracts. High fidelity gates with small step size and a high angle were necessary to accurately capture the highly curved anatomy of the SM.

#### Gate Placement

Following consultation with the Department of Anatomy at Trinity College Dublin, a Boolean logic protocol was designed to isolate the SM using AND and OR gates within ExploreDTI. Emphasis was given towards gates that: (a) included the maximum possible number of SM streamlines; (b) were positioned around landmarks clearly identifiable on DW or T1 images; and (c) that would consistently capture at least some section of the SM. Initial gate consensus consisted of two OR gates in the axial plane at the agreed anticipated anterior and posterior margins of the SM and an AND gate in the coronal plane as the SM arcs over the inter-thalamic adhesion. Three AND gates would generate the same final result, however in our exploratory investigations and rater training session, two OR gates were selected at the anterior and posterior margins, so our protocol continued to use this for consistency.

As the SM neurons merge into a solitary tract just behind the anterior commissure (**Figure 1**), the uppermost section of the commissure as seen on DWI was used as the primary landmark for placement of the anterior axial OR gate (**Figure 2A**). The SM terminates posteriorly in both habenulae (**Figure 1**). This tiny structure is not seen consistently in either DW or T1 images (Lawson et al., 2013). The habenulae bilaterally are within the stalks of the pineal gland, a structure consistently seen on T1 MRI. However, the position of the pineal relative to its stalk as seen on MRI depends on the size of the gland and on how the pineal ''hangs'' within the quadrigeminal cistern. To reliably include part of the SM tract, the superior most point of the pineal on T1 imaging was established as the only consistent primary landmark for placement of the posterior OR gate (**Figures 1**, **2B**). Occasionally this resulted in a shortened SM that omitted the posterior end. Finally, the SM always arcs between both thalami and above the inter-thalamic adhesion if present (**Figure 1**). The inter-thalamic adhesion is not present in 20%–30% of individuals (Allen and Gorski, 1991); as such the inner half of both thalami on coronal T1 imaging was used as the primary landmark for the final AND gate (**Figure 2C**). Gate placement is described in detail in **Figure 2**.

Following an initial training period on a different dataset of 10 subjects, two raters placed the gates on DW and T1 MR images of 50 subjects. Raters were blind to age and sex of the subjects. The gates were applied to the low fidelity scout whole brain tracts to generate streamlines of interest encompassing the SM. Larger obviously non-SM tracts (most notably the fornix) were excluded using NOT gates; see **Figure 4** for details of typical NOT gates. The remaining streamlines were split and segmented (using the Splitter and Segmentation tools in ExploreDTI) to include all fibers that only travel between the three gates yet eliminate those sections of the streamlines peripheral to the anterior and

FIGURE 6 | Superior view of a bilateral SM on an axial T1 plane. Note its posterior insertion into the high signal intensity regions at the posterior thalamus, consistent with white matter. Male, aged 27.

posterior gates (see **Figure 5**). Once satisfied that the SM could be visualized within the section, the high fidelity whole brain tracts were reloaded with the gates in place. The low fidelity scout tracts were then discarded. With the gates now in place, the high fidelity whole brain tracts could be used to ''dissect'' out the region around the SM.

#### Cleanup

The next step required that streamlines not associated with the SM be removed. Cleanup was achieved with each rater and a neuroanatomist examining every subject together, collectively excluding streamlines not considered SM fibers using NOT gates. Criteria for exclusion included streamlines deviating from


The effect of gender on each measure following ANCOVA, correcting for age and eTIV is shown. The effect of age on each measure following partial correlation analysis correcting for sex and eTIV is also shown. AD, Axial Diffusivity; ANCOVA, Analysis of Covariance; CL, Linear Diffusion; CI, Confident Interval; CP, Planar Diffusion; CS, Spherical Diffusion; FA, Fractional Anisotropy; ICC, Inter-rater correlation coefficient; MD, Mean Diffusivity; RD, Radial Diffusivity.

the known arch shape of the SM and streamlines consistent with known adjacent tracts. To address potential bias in the same subject between raters, the neuroanatomist was blind to the identity of each subject. These NOT gates were variable according to each subject and depended on the individual neuroanatomy of the region. However, commonly excluded streamlines included fibers from the body of the fornix, the stria terminalis, the median forebrain bundle, superior thalamic peduncle, and various inter-thalamic connections. A second neuroanatomist reviewed a random subsample of 25 tracts from both raters to check face validity of the SM (looking at start and end points, the arch shape, the tract's relationship to other structures on MRI images (thalamus, fornix, anterior/posterior commissures and third ventricle) as well as the overall general impression of the tract. The second neuroanatomist also checked both raters' tracts simultaneously to investigate if they were occupying approximately the same space. This was considered a sufficient subsample for quality control.

Tract statistics including dimensions (length and volume), standard diffusion measures such as FA, mean diffusivity (MD), axial diffusivity (AD) and RD and Westin diffusion measures [Spherical Diffusion (CS), Linear Diffusion (CL) and Planar Diffusion (CP)] were calculated on each tract using ExploreDTI's built in tract statistics plugin.

### Estimated Total Intracranial Volume and Statistical Analysis

Estimated Total Intracranial Volume (eTIV) was calculated from the T1 images using the Freesurfer 6.0 image analysis suite<sup>3</sup> (Fischl, 2012). The technical details of this procedure are described elsewhere (Buckner et al., 2004; Voevodskaya et al., 2014; Sargolzaei et al., 2015).

<sup>3</sup>http://surfer.nmr.mgh.harvard.edu/

FIGURE 8 | Close up side view of a single right-sided SM on a T1 sagittal plane (the left side is hidden behind the T1 sagittal plane) showing the stria arching over the inter-thalamic adhesion. Male, aged 40.

#### Statistics

A comparison of reliability between each independent rater was performed using interclass correlation coefficient (ICC) analysis using SPSS24<sup>4</sup> for the following measurements: dimensional measures, SM length and tract volume, standard diffusion metrics, FA, MD, AD and RD and Westin diffusion measures, CS, C<sup>L</sup> and CP. Rater measures were averaged to represent a mean DWI measure for each metric and the effect of gender on each diffusion metric was determined using independent Analysis of Covariance (ANCOVA) in SPSS correcting for age and eTIV. Finally, partial correlation analyses were performed to examine the effect of age on each diffusion metric when controlling for gender and eTIV.

#### RESULTS

The SM was generated in 50 subjects by two raters and reconstructed consistently in every subject (**Table 1**). Even though the SM can be considered a bilateral structure, both halves unite into a single tract as it makes its way from front to back. As such it was not possible to reliably separate out the left from right SM using tractography. These results are presented as a single total bilateral (i.e., left and right SM combined) tract. Examples of generated SM are found in **Figures 6**–**10**. **Figure 10** demonstrates some of the inter-subject variability seen in the tract. A video showing a rotating SM to demonstrate the structure in 3D is included in the supplementary section. The second neuroanatomist had concerns with one generated SM (the arch appeared flattened on inspection), which was redone afresh and cleaned as per protocol. Both neuroanatomists and rater agreed with the changes.

As expected, the structure arched dorso-caudally until about mid-thalamus and then continued ventro-caudally as it neared the habenula. Where there was an inter-thalamic adhesion, the SM arched dorsal to this. A generated tract in relation to anatomical landmarks is shown in **Figure 11**. Our current protocol captures the SM only from the anterior commissure rostrally to where the arch reaches the level of the uppermost point of the pineal gland as it continues to arch ventrally (see **Figure 1**, and for more details see ''Materials and Methods'' section). This landmark was chosen for anatomical and imaging consistency, even though it routinely resulted in a slightly shorted SM, omitting the caudal most aspect of the tract where it meets the habenula. The mean lengths of the tract were 29.8 mm and 30.3 mm for Rater 1 and 2, respectively. The mean volumes of the SM were 129 mm<sup>3</sup> and 131 mm<sup>3</sup> for Rater 1 and 2. Inter-subject variability is shown in **Table 2**. There is no reported SM length in the literature to our knowledge. However, these lengths correspond with measured SM lengths in eight cadaveric brains by the Anatomy Department at Trinity College Dublin of between 28 mm and 35 mm.

Good inter-rater reliability was demonstrated for each measurement (mean ICC 0.91 over all measures). There were no significant differences between males and females for any DWI metric examined. Partial correlation analyses between diffusion metrics and age identified a positive correlation for RD only

<sup>4</sup>https://www.ibm.com/analytics/us/en/technology/spss/

FIGURE 9 | An example of a left and right fused SM on axial and coronal T1 image planes. The left and right sides are fused in the middle of the tract. Note its posterior insertion into the high signal intensity regions at the posterior thalamus, consistent with white matter. Female, aged 22.

(p = 0.038, r = 0.3) when correcting for gender and eTIV. No other metric showed a significant correlation with age (p < 0.05); refer to **Tables 1**, **2** for results.

#### DISCUSSION

This study describes a reliable technique for reconstructing the functionally important SM using diffusion-weighted tractography. The technique demonstrates excellent dimensional inter-rater reliability between two independent raters and face validity from a second independent neuroanatomist. Tract metrics also showed excellent inter-rater reliability with high ICCs for dimensional measures (length and volume) and diffusion metrics. The mean tract length were 29.8 mm (Rater 1) and 30.3 mm (Rater 2) and correspond to observed tract lengths of between 28 mm and 35 mm in cadaveric brains. Of note, there was no difference in SM dimensions or diffusion metrics between males and females.

Only one SM diffusion metric changed with age. RD, a measure of diffusivity perpendicular to the fiber tract orientation increased with age. AD, a measure of diffusivity parallel to the tract, showed no change with age. FA showed no change with age. This is a measure of the relative diffusivity along the bundle. Often in aging white matter tracts, FA decreases as RD increases (Bennett et al., 2010). This is due to the relatively greater change of diffusivity along the tract compared to across it. This is not the case in our SM reconstructions, where only RD increases. The increase in RD with age is consistent with other studies reporting an increase across multiple white matter brain regions (Hsu et al., 2008) and particularly in limbic areas (Kumar et al., 2013). This RD change may reflect deteriorations in myelin with age (Song et al., 2002; Wu et al., 2011), and myelin loss is considered a normal part of healthy aging (Peters, 2002). Dorsal diencephalic conduction system input changes through an aging SM may manifest as overall changes in mood, hedonic states and motivation that occur with age (Prince et al., 1999; Harada et al., 2013; Sutin et al., 2013).

Only one previous study has investigated the SM using diffusion-weighted imaging (Kochanski et al., 2016). This study used a probabilistic local approach in five pre-surgical patients defining the estimated site of the habenula as a sole seed

point. Using a single probabilistic seed point with a local approach may have limitations in capturing the known gross anatomy of the SM. Indeed the tracts in this study do not arch over the inter-thalamic adhesion, unlike our SM tracts. A local seed based approach is susceptible to errors as the tract propagates (Reisert et al., 2011; Anastasopoulos et al., 2014). We used three anatomically guided gates to accurately identify the known arch-like shape, followed by careful cleaning according to known boundary structures. This anatomically guided approach on global (whole brain) tracts allowed us to predefine our tract start (anterior) and end (posterior) points for consistency. Having tracts with defined boundaries will allow easy comparison between individuals for future studies. Using a third gate across the thalamus and above the interthalamic adhesion allowed us to focus solely on streamlines that corresponded to the known anatomical arch of the SM and reduced the likelihood of generating deviant tracts. We believe our approach more accurately describes the trajectory of the SM and could be useful in conjunction with standard structural imaging in perioperative stereotactic localization of the SM for DBS. Furthermore, the habenula is a small structure of approximately 30–36 mm<sup>3</sup> (Savitz et al., 2011; Lawson et al., 2013) and may not be accurately identified on standard MRI especially at 3T (Strotmann et al., 2014). We felt it more appropriate to use an imaging fixed point at the pineal to optimize consistency and reliability. Our present study has some other advantages including a higher number of subjects (50), healthy controls instead of patients, prior neuroanatomical planning and face validity checks from a second neuroanatomist.

#### Strengths and Limitations

Deterministic tractography is an experimental method that is subject to some technical challenges (Jones, 2008). There are known issues involving curving and kissing tracts. The SM is highly curved and turns 180◦ in a relatively tight space. This was addressed by reconstructing the SM using high angle (89◦ ), small step-size (0.5 mm) whole brain tracts. The SM lies deep within the center of the brain, surrounded by and touching other white matter bundles; the fornix, stria terminalis and various thalamic connections. CSD was chosen as the most appropriate tractography technique due to its superior ability at detecting complex crossing and kissing fibers in a recent analysis of fiber estimation approaches (Wilkins et al., 2015). As with all diffusion-weighted reconstruction, the generated SM only describes the diffusivity in the region of interest and it is assumed that this corresponds to an actual tract. The length of our generated SM (30 mm) approximates to the known length of the SM. Also, a subsample of images (50%) was presented to a second neuroanatomist to check face validity of the results. To address this further we plan to systematically measure the length and curvature of a number of cadaveric SMs to determine baseline dimensions of this tract for comparison.

An experienced neuroanatomist with expertise of the region aided the rater during cleanup. We felt this represented the

most valid approach for removing extraneous streamlines. As the protocol was designed to encompass as much of the SM as possible in a consistent manner, the initial gates were overly broad. This resulted in a large amount of cleanup. The average time of cleanup (both rater and neuroanatomist working together) was estimated at about 20 min per subject. It is appreciated that this approach may not be available or practical for some centers. However, with adequate neuroanatomical training, this limitation may be overcome and a researcher/technician experienced with ExploreDTI and the complex regional anatomy could single-handedly undertake the entire process.

One limitation is that the posterior limit of the tract was occasionally clipped due to the to position of the gate at the pineal gland. The terminus of the SM anatomically is the habenula, which lies within the stalks of the pineal. The habenulae/pineal stalks are inconsistently seen on standard T1 images (Lawson et al., 2013), and as such the pineal was chosen as a marker. However, we found during our initial exploration of potential gates that only a gate placed just above the level of the most superior point of the pineal guaranteed inclusion of the SM. Our reconstructed tracts are therefore slightly shorter than anatomical ones. Higher resolution T1 imaging showing the pineal stalk may address this problem for future studies.

#### TABLE 2 | Inter-subject SM dimension variability.


Volumes, lengths and standard errors are shown. The upper and lower limits of data within the 95% confidence interval are also shown.

The mean FA (0.251) calculated for our tract was low compared to other limbic tracts e.g., stria terminalis (FA = 0.37; Kwon et al., 2011), mammillothalamic tract (FA = 0.38; Kwon et al., 2010), pre-commissural (FA = 0.341) and post commissural fornix (FA = 0.359; Yeo et al., 2013). This could be explained by the fact the SM exists as a small solitary tract within the third ventricle surrounded by cerebrospinal fluid (CSF) on three sides for most of its course. The diameter of the SM as measured through our initial cadaveric measurements was between 1.5 mm and 2.5 mm. As our voxel size was 2 mm<sup>3</sup> , it is probable that partial volume effects where both CSF and SM mixed within voxels were present in this region. Partial voluming from inclusion of CSF has been shown to collapse FA (Alexander et al., 2001). FA values of below 0.25 in otherwise healthy adults with normal appearing white matter may reflect this effect rather than changes in white matter coherence (Pfefferbaum and Sullivan, 2003). We plan to address this in future research by decreasing our voxel size to 1 mm<sup>3</sup> and also through using free water correction methods (Pasternak et al., 2009).

#### Implications

As the first component of the dorsal diencephalic conduction system, the SM is fundamental in the flow of information from frontal-limbic regions to key midbrain monoaminergic nuclei (Sutherland, 1982). The habenula has been shown to be involved in aversive conditioning, reward pathways (Boulos et al., 2017) and in the pathophysiology of several psychiatric conditions (Fakhoury, 2017). The DCCS seems to be particularly relevant in negative reward prediction errors, or ''disappointment'', where failure to receive an expected reward integrates cognitive and emotional inputs through the habenula to affect a depressive behavioral response (Kaye and Ross, 2017).

In particular, changes in structure (Savitz et al., 2011; Schmidt et al., 2017) and function (Furman and Gotlib, 2016; Lawson et al., 2017) of the habenula have been shown in depression. This makes sense considering the DDCS's pivotal role in mediating serotonergic tone. Depressive symptoms, in particular anhedonia and motivation, have been associated with disrupted habenula function (Lawson et al., 2017; Liu et al., 2017). This reflects its role in connecting the basal septal pleasure areas and midbrain monoaminergic response. However, the habenula's main afferent input from these frontal regions, the SM, has not been a focus of clinical study, mainly due to the difficulties in assessing the tract in vivo. Using our approach, it may be possible to use diffusion metrics to investigate the SM in a wide range of neuropsychiatric conditions such as depression and anxiety (serotonin), schizophrenia and addictions (dopamine) and pain (noradrenaline). Changes in afferent habenula inputs may be reflected in diffusion metric changes, focusing attention on regions upstream to the DCCS such as the septal pleasure areas. Future research using higher resolution imaging may also reveal the lateral and medial habenular components (Strotmann et al., 2014) and their discrete afferents via the SM using diffusion connectomics, furthering our understanding of the complex and distinct connections between the forebrain and midbrain.

The SM has long been suggested as a potential DBS target in treatment-resistant depression (Blander and Wise, 1989; Sartorius and Henn, 2007; Sartorius et al., 2010; Kiening and Sartorius, 2013; Kochanski et al., 2016). This is due its key role in integrating diverse basal forebrain into a solitary white matter bundle prior to processing in the habenula. It has been attempted in one patient who showed remarkable improvement from intractable depression following high frequency stimulation of the SM (Sartorius et al., 2010). Our tractography approach may aid future 3D stereotactic localization of the SM prior to electrode placement in such a patient, as the habenulae and SM are inconsistently seen on standard structural MRI (Lawson et al., 2013; Kochanski et al., 2016). As well as providing another layer of anatomical certainty in combination with standard imaging during electrode insertion, SM diffusion metrics may also offer insights into the microstructural integrity of the tract pre and post treatment.

This study establishes a novel and reliable method for reconstructing the SM using deterministic diffusion-weighted tractography The SM was reconstructed using high angle, small step-size tracts to capture the particularly small and highly angular path of the SM. The use of neuroanatomical expertise during the planning, cleanup and checking processes strengthened the face validity of the technique. Reconstructing the SM in a reliable and accurate fashion may contribute to research into neuropsychiatric disorders and treatments.

### AUTHOR CONTRIBUTIONS

DWR: recruitment, protocol design, data analysis, manuscript preparation and neuroanatomy. ER: recruitment, protocol design, tractograpy, data analysis and manuscript preparation. SR and SA: recruitment and tractography. CF, KD and LT: recruitment. KJL: data analysis and neuroanatomy. DB and PT: protocol design and neuroanatomy. TF: principal investigator. VO: principal investigator and manuscript preparation. EO: recruitment, protocol design, data analysis and manuscript preparation.

### FUNDING

This project was funded by the Irish Health Research Board as part of the REDEEM (Research in Depression, Endocrinology, Epigenetics and neuroiMaging) study at Trinity College Institute of Neuroscience and the Department of Psychiatry, Trinity College Dublin (TCD). Grant code: 201651.12553.

### ACKNOWLEDGMENTS

We acknowledge the assistance of Trinity College Institute of Neuroscience especially Mr. Sojo Joseph. We also acknowledge the individuals who donated their bodies to the anatomy department at TCD, without who's help we would not have been able to develop this study.

#### REFERENCES


Yeo, S. S., Seo, J. P., Kwon, Y. H., and Jang, S. H. (2013). Precommissural fornix in the human brain: a diffusion tensor tractography study. Yonsei Med. J. 54, 315–320. doi: 10.3349/ymj.2013.54.2.315

**Conflict of Interest Statement**: 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.

Copyright © 2018 Roddy, Roman, Rooney, Andrews, Farrell, Doolin, Levins, Tozzi, Tierney, Barry, Frodl, O'Keane and O'Hanlon. 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 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.

# A Quantitative Tractography Study Into the Connectivity, Segmentation and Laterality of the Human Inferior Longitudinal Fasciculus

Sandip S. Panesar<sup>1</sup> , Fang-Cheng Yeh1,2, Timothée Jacquesson<sup>1</sup> , William Hula<sup>3</sup> and Juan C. Fernandez-Miranda<sup>1</sup> \*

<sup>1</sup> Department of Neurological Surgery, University of Pittsburgh, Pittsburgh, PA, United States, <sup>2</sup> Department of Bioengineering, University of Pittsburgh, Pittsburgh, PA, United States, <sup>3</sup> Veterans Affairs Pittsburgh Healthcare System, Pittsburgh, PA, United States

The human inferior longitudinal fasciculus (ILF) is a ventral, temporo-occipital association tract. Though described in early neuroanatomical works, its existence was later questioned. Application of in vivo tractography to the neuroanatomical study of the ILF has generally confirmed its existence, however, consensus is lacking regarding its subdivision, laterality and connectivity. Further, there is a paucity of detailed neuroanatomic data pertaining to the exact anatomy of the ILF. Generalized Q-Sampling imaging (GQI) is a non-tensor tractographic modality permitting high resolution imaging of white-matter structures. As it is a non-tensor modality, it permits visualization of crossing fibers and accurate delineation of close-proximity fiber-systems. We applied deterministic GQI tractography to data from 30 healthy subjects and a large-volume, averaged diffusion atlas, to delineate ILF anatomy. Post-mortem white matter dissection was also carried out in three cadaveric specimens for further validation. The ILF was found in all 60 hemispheres. At its occipital extremity, ILF fascicles demonstrated a bifurcated, ventral-dorsal morphological termination pattern, which we used to further subdivide the bundle for detailed analysis. These divisions were consistent across the subject set and within the atlas. We applied quantitative techniques to study connectivity strength of the ILF at its anterior and posterior extremities. Overall, both morphological divisions, and the un-separated ILF, demonstrated strong leftwardlateralized connectivity patterns. Leftward-lateralization was also found for ILF volumes across the subject set. Due to connective and volumetric leftward-dominance and ventral location, we postulate the ILFs role in the semantic system. Further, our results are in agreement with functional and lesion-based postulations pertaining to the ILFs role in facial recognition.

Keywords: inferior longitudinal fasciculus, tractography, semantic anatomy, white matter anatomy, non-tensor tractography

### INTRODUCTION

The human inferior longitudinal fasciculus (ILF) is a ventral, temporo-occipital association tract. Historically, it was described as connecting the superior, middle, inferior and fusiform gyri, to the lingual, cuneate, lateral-occipital and occipito-polar cortices. Early descriptions came from post-mortem white matter dissection (Dejerine and Dejerine-Klumpke, 1895; Crosby, 1962;

#### Edited by:

Laurent Petit, Centre National de la Recherche Scientifique (CNRS), France

#### Reviewed by:

Francesco Latini, Uppsala University, Sweden Alessandro De Benedictis, Ospedale Pediatrico Bambino Gesù (IRCCS), Italy

#### \*Correspondence:

Juan C. Fernandez-Miranda fernandezmirandajc@upmc.edu

> Received: 24 March 2018 Accepted: 18 May 2018 Published: 05 June 2018

#### Citation:

Panesar SS, Yeh F-C, Jacquesson T, Hula W and Fernandez-Miranda JC (2018) A Quantitative Tractography Study Into the Connectivity, Segmentation and Laterality of the Human Inferior Longitudinal Fasciculus. Front. Neuroanat. 12:47. doi: 10.3389/fnana.2018.00047

**63**

Gloor, 1997). More recently, radioisotopic tracer studies in non-human primates questioned the existence of a robust temporo-occipital fasciculus, proposing instead a series of cortico-cortical U-fibers, traveling along the antero-posterior temporal distance (Tusa and Ungerleider, 1985; Catani et al., 2003). The introduction of tractography somewhat reconciled this controversy: Early diffusion tensor imaging (DTI) studies confirmed the existence of the human ILF (Catani et al., 2002, 2003; Mori et al., 2002; Jellison et al., 2004; Kier et al., 2004; Wakana et al., 2004; Fernández-Miranda et al., 2008; Urbanski et al., 2008). Catani et al. (2002, 2003) described a bicomponent structure, with a direct temporo-occipital sub-tract, and an indirect series of cortico-cortical U-fibers. The direct subfascicle demonstrated a tri-pronged posterior, and bi-pronged anterior connectivity profile. This connective and structural description remained largely unchanged until a recent dissection and DTI series by Latini (2015) and Latini et al. (2017). The authors proposed a rightward-lateralized ILF arrangement derived from its posterior connectivity profile: dorsally, distinct ILF subfasciculi originated from the cuneate and lateral occipital lobes, while ventrally, subfasciculi originated from the fusiform and lingual gyri (Latini, 2015; Latini et al., 2017).

Prior to the introduction of functional neuroimaging, data pertaining to ILF function was derived primarily from lesional data. Visual agnosia and amnesia were attributed to bilateral damage to ventral temporal white matter. Prosopagnosia (Benson, 1974; Bauer, 1984; Michel et al., 1989; Grossi et al., 2014) and visual hypo-emotionality (Habib, 1986; Takeuchi et al., 2013; Fischer et al., 2016) were linked to right-sided ventral temporal lesions. Further insights have come from both functional neuroimaging (e.g., fMRI) and intraoperative electrical stimulation (IES) (Duffau et al., 2005; Mandonnet et al., 2007). fMRI studies demonstrated semantic activations within temporo-polar, anterior temporal, basal occipito-temporal and occipital lobes corresponding to dominant-hemisphere ILF trajectory (Saur et al., 2008, 2010). An IES study by Mandonnet et al. (2007) proposed the ILF's role in semantic functionality, consistent with the 'dorsal-ventral' model of language organization (Hickok and Poeppel, 2015).

Tractography permits functional insights by studying corticalconnectivity patterns of white matter tracts. It also allows calculation of volumetric, lateralization and subdivision patterns of white matter systems. It therefore surpasses post-mortem dissection in terms of ability and accuracy (Fernandez-Miranda et al., 2012). DTI is unable to delineate crossing fibers and demonstrate cortical connectivity accurately (Fernandez-Miranda et al., 2012; Farquharson et al., 2013). Generalized q-sampling imaging (GQI) dispenses with the diffusion tensor, allowing crossing fibers at close proximity to be tracked at highresolution (Yeh et al., 2010, 2013). Our group has pioneered the application of GQI tractography to neuroanatomical (Wang et al., 2013, 2016; Fernández-Miranda et al., 2015; Meola et al., 2016; Yoshino et al., 2016; Panesar et al., 2017) and surgical studies (Fernandez-Miranda et al., 2012). Further, our use of large-volume GQI-derived atlases (Yeh and Tseng, 2011; Wang et al., 2013; Fernández-Miranda et al., 2015; Panesar et al., 2017) provides a model of 'average' white matter anatomy for validation. With these considerations in mind, we set out to study the subdivision, asymmetry and connectivity of the human ILF using GQI tractography using quantitative connectometry methods, along with diffusion atlas and dissection validation.

## MATERIALS AND METHODS

### Participants

We conducted a subject-specific deterministic fiber tractography study in 30 right-handed, neurologically healthy male and female subjects, age range 23–35. The data were from the Human Connectome Project (HCP) online database [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. Likewise, data from 1021 individual HCP subjects were utilized to compile the averaged diffusion atlas.

### Image Acquisition and Reconstruction

The HCP diffusion data for individual subjects were acquired using a Siemens 3T Skyra system, with a 32-channel head coil (Siemans Medical, Erlangen, Germany). A multi-shell diffusion scheme was used, and the b-values were 1000, 2000, and 3000 s/mm<sup>2</sup> . The number of diffusion sampling directions were 90, 90, and 90, respectively. The in-plane resolution and slice thickness were both 1.25 mm (TR = 5500 ms, TE = 89 ms, resolution = 1.25 mm × 1.25 mm, FoV = 210 mm × 180 mm, matrix = 144 × 168). The DSI data were reconstructed using the generalized q-sampling imaging approach (Yeh et al., 2010) using a diffusion distance ratio of 1.2 as recommended by the original study.

### HCP 1021 Atlas

A total of 1021 participants from the HCP database were used to construct the atlas. The image acquisition parameters are identical to those described previously. The diffusion data were reconstructed and warped to the Montreal Neurological Institute (MNI) space using q-space diffeomorphic reconstruction (Yeh and Tseng, 2011) with a diffusion sampling length ratio of 1.25 and output resolution of 1 mm. The group average atlas was then constructed by averaging the reconstructed data of the 1021 individual subjects within the MNI space.

## Fiber Tracking and Analysis

We performed deterministic fiber tracking using DSI Studio software, which utilizes a generalized streamline fiber tracking method (Yeh et al., 2013). Parameters selected for fiber tracking included a step size of 0.2 mm, a minimum fiber length of 20 mm and a turning angle threshold of 60◦ . For progression locations containing >1 fiber orientation, fiber orientation most congruent with the incoming direction and turning angle <60◦ was selected to determine subsequent moving direction. Each progressive voxels' moving directional estimate was weighted by 20% of the previous voxels incoming direction and by 80% of its

nearest fiber orientation. This sequence was repeated to create fiber tracts. Termination of the tracking algorithm occurred when the quantitative anisotropy (QA) (Yeh et al., 2013) dropped below a subject-specific value: when fiber tract continuity no longer met the progression criteria, or when 100,000 tracts were generated. We pre-selected QA termination threshold, between 0.02 and 0.08, by analyzing the number of false continuities generated within each subjects' dataset and chose the compromise value that allowed optimal anatomical detail with minimal noise. Likewise, we selected a smoothing parameter of 50% for the same reason stated previously.

### Region of Interest Placement and Fiber Selection

Unlike our previous studies, and in an effort to replicate published methodology (Catani et al., 2003; Latini et al., 2017), we chose a two region of interest (ROI) approach rather than an atlas-based seeding approach. This method was chosen to minimize a priori selection of feed-backward or feed-forward tracts. A spherical ROI was positioned within the white matter of the anterior temporal lobe. A second, rectangular ROI was placed in the coronal plane at the temporo-occipital junction. To avoid commissural fibers and false continuities, a rectangular region of avoidance was placed across the sagittal hemispheric midline. Once 100,000 fiber tracts were generated, we manually removed fibers passing through the ventral external capsule [i.e., inferior fronto-occipital (IFOF), uncinate (UF) fasciculi] and fibers resembling the arcuate fasciculus (AF), superior longitudinal fasciculus (SLF) and middle longitudinal fasciculus (MdLF) (**Figure 1A**).

#### Quantitative Connectivity Analysis and Connectogram Creation

We used a quantitative method to define cortical terminations. Each sub-fascicle was merged with its hemispheric counterpart (i.e., left and right). The 'connectivity matrix' function in DSI studio was used to generate matrices representing the number of fibers terminating within regions of a specially modified, per-subject aligned version of the Automated Anatomical Labeling (AAL) atlas. Three connectivity matrices were generated per subject, to give a total of 90 connectivity matrices corresponding to whole ILF (i.e., undivided), dorsal ILF and ventral ILF (see section "Results") over the 30 individual subjects. These matrices were collated in Matlab and the number of fibers corresponding to each respective connection were divided by the total number of fibers for each fascicle, per subject. Values were then scaled over the range of 30 subjects to give a connection index (CI) between 0 and 100, with 0 representing no connectivity, and 100 representing strongest relative connectivity to a particular atlas region. Index values were then used to generate connectograms for each sub-fascicle and for whole ILF bundles, using CIRCOS<sup>1</sup> . Connectograms provide a unique method of visualizing network topology by demonstrating weighted strength of connectivity between brain regions.

### Volumetry and Lateralization

We calculated the number of voxels occupied by each fiber trajectory (streamlines) and the subsequent volume (in milliliters) of each gross ILF bundle and their manually separated subfascicles. The volumes of left and right merged ILFs and sub-fascicles across 30-subjects were each subjected to an independent samples T-test in SPSS (IBM Corporation, Armonk, NY, United States) to calculate significance of mean hemispheric volumetry over the 30 subjects. We also calculated the lateralization index (LI) using the formula {[(volume left − volume right)÷(volume left + volume right)] × 2] (Catani et al., 2007).

#### White Matter Dissection

Three human hemispheres were prepared for dissection. First, they were fixed in a 10% formalin solution for 2 months. After fixation, the arachnoid and superficial vessels were removed. The brains were subsequently frozen at −16◦C for 2 weeks, as per the Klingler method (Ludwig and Klingler, 1956). Dissection commenced 24 h after the specimens were thawed and proceeded in a step-wise, superficial-to-deep process, at the lateral and inferior surface of the temporal lobe. Dissection was achieved using wooden spatulas to remove successive layers of gray matter followed by white matter. We utilized post-mortem white matter dissection only to study the ILFs' position amongst other large association and projection white matter tracts, rather than for connectivity reasons.

### RESULTS

### Morphology and Subdivision of ILF Bundles

Bundles resembling the ILF were found in 60/60 hemispheres (**Figure 1B**) in our subject-specific tractographic analysis and bilaterally within the 1021-subject template. Upon generating and extracting whole ILF bundles we observed a distinct posterior termination pattern across the individual subjects and in the template, which we used to further divide the ILFs (**Figure 1C**). Though the whole ILF bundles demonstrated a distinct posterior termination morphology, their anterior aspects did not. Upon morphological assessment, it was clear that the ILF followed a distinctive 'dorsal-ventral' termination pattern. We thus chose to manually separate these components, terming them 'dorsal' and 'ventral,' in addition to assessing the whole, unseparated bundle (**Figure 1D**). Both dorsal and ventral bundles were found in 60/60 of individual subjects' hemispheres. Left and right ILFs were successfully replicated on the 1021-subject atlas and assumed identical morphology and subdivisions compared to those from the subject-specific study (**Figures 2A–D**).

#### Quantitative ILF Connectivity Dorsal ILF Sub-Fascicle (Table 1)

The left dorsal ILF bundle demonstrated strongest connectivity patterns between the superior occipital to superior temporal gyri (CI: 100.00) and middle temporal gyri (CI: 41.13). All other

<sup>1</sup>http://mkweb.bcgsc.ca/tableviewer/visualize/

connections had a CI <40. Right sided dorsal ILF connectivity was comparatively weaker versus the left: maximum connectivity was found between the cuneus and superior temporal gyrus (CI: 25.99), with all other connections demonstrating CI <25. Left dorsal ILF connectivity was spread over nine occipito-temporal connections, while the right was spread over 7. The average CI for the left dorsal ILF was 14.77 (SD = 26.13), while it was 10.45 (SD = 15.83) for the right (Not significant, T = 0.61, p = 0.547). Median CI for the left dorsal ILF was 23.10 and was 5.90 on the right (**Figure 3**).

#### Ventral ILF Sub-Fascicles (Table 2)

The left ventral ILF demonstrated strongest connectivity between the lingual gyrus and middle temporal gyrus (CI: 100.00), followed by lingual to superior temporal gyrus (CI: 74.63), and lingual to inferior temporal gyrus (CI: 64.86). Left ventral ILF also demonstrated strong connectivity between the calcarine area to middle temporal (CI: 63.17), superior temporal (CI: 58.29) and inferior temporal (CI: 55.24) gyri. All other left-sided ventral connections had CI's under 50. The right ventral ILF demonstrated relatively weaker connectivity than the left, but with similar patterns: Lingual gyrus to middle temporal (CI: 88.66), inferior temporal (CI: 48.77) and superior temporal (CI: 41.34) gyri. All other connections had CI less than 40. Left ventral ILF demonstrated 22 distinct connections, whilst the right demonstrated 21. The average CI for the left ventral ILF was 16.76 (SD = 25.80) and 8.93 (SD = 17.10) for the right (not significant, T = 1.64, p = 0.053). Median CI for the left ventral ILF was 3.29. It was 0.78 for the right (**Figure 4**).

#### Whole ILF Bundles (Table 3)

When taken as whole, merged bundles, both left and right ILFs each demonstrated 22 patterns of occipito-temporal connectivity, respectively. The left ILF, however, demonstrated stronger connectivity indices versus the right. On the left, strongest connectivity was demonstrated between lingual to middle temporal gyrus (CI: 100.00), superior occipital to superior temporal gyrus (CI: 77.81), lingual to superior temporal gyrus (CI: 77.67), calcarine to middle temporal gyrus (CI: 69.94) and calcarine to superior temporal gyri (CI: 63.18). All other connections had CI <60. The strongest connections on the right side were between lingual and middle (CI: 77.54), inferior (CI: 41.05) and superior (CI: 40.70) temporal gyri. All other connections had a CI <40. The average CI for the whole ILF was 29.50 (SD = 29.29) on the left, while it was 10.60 (SD = 15.83) for the right (significant, T = 3.07, p = 0.004). Median CI for the left ILF was 23.10 and was 5.90 for the right (**Figures 5**, **6A**).

#### Volumetry and Lateralization

Whole, merged ILF bundles had a mean volume of 19.1 ml (SD = 6.0), while right sided whole ILFs had a mean volume of 14.1 ml (SD = 5.4). This 5 ml difference was significant (T = 3.38,


Darker color represents stronger index.

p = 0.001). For separated sub-fascicles, mean left dorsal ILF volume was significantly larger than right; at 9.8 ml vs. 6.9 ml, respectively (T = 3.50, p = 0.001). For the ventral ILF, mean left sided volume was 14.7 ml (SD = 5.5). For the right ventral ILF it was 10.9 ml (SD = 5.6). This hemispheric difference was significant (T = 2.64, p = 0.005). For whole ILF bundles, LI was 0.31 demonstrating leftward asymmetry. Both dorsal and ventral ILFs' had a mean LI of 0.3, respectively. Using the HCP 1021 atlas, the whole ILF was 17.7 ml on the left and 11.0 ml on the right. When divided, the left dorsal ILF was 8.8 ml and the left ventral ILF was 11.1 ml. The right dorsal ILF was 2.4 ml, whilst the right ventral ILF was 9.4 ml.

### White Matter Tract Relations as Revealed by Tractography

At its anterior aspect, The ILF shares terminations within the temporal pole with the ventral UF. It then travels posteriorly, lateral to the temporal horn. At its longitudinal-middle segment, temporal connections of the AF overlie it. Originating from the superior temporal gyri, the MdLF passes deep to the posterior limit of the Sylvian fissure and the AF, traveling obliquely to the dorsal occipital and superior parietal area, remaining dorsal to the ILF at all stages. At the temporo-occipital junction, the ILF travels posteriorly, joining the sagittal stratum inferiorly and remaining ventro-lateral to the IFOF. As such, the IFOF lies between the ILF and optic radiations, all of which traverse anteroposteriorly via the sagittal stratum to occipital terminations. As it passes over the basal temporal areas, i.e., the fusiform gyrus, the ILF is overlaid by the vertical occipital fasciculus (VOF), a short vertically oriented white matter tract originating from the basal occipito-temporal areas and connecting with the dorsolateral occipital gyri. At its posterior and occipital aspect, the ILF shares termination areas with the IFOF, MdLF, VOF and the optic radiations (**Figures 6B,C**, **7A,B**).

### White Matter Relations as Revealed by Dissection

From the lateral surface, the dissection started at the intersecting region of the parieto-occipito-temporal junctions, removing the U-fibers overlying the supramarginal and angular gyri. The

FIGURE 3 | A connectogram representing bilateral connectivity patterns of the dorsal ILF subfascicle. Atlas regions are listed around the circumference of the connectogram. TInf, inferior temporal gyrus; TMd, middle temporal gyrus; TSp, superior temporal gyrus; Fu, fusiform gyrus; OInf, inferior occipital gyrus; OMd, middle occipital gyrus; OSp, superior occipital gyrus; Cal, calcarine gyrus; Cun, cuneus; Lin, lingual gyrus. Left and Right hemispheric connections are demonstrated with a suffix \_L or \_R.



Darker color represents stronger index.

FIGURE 4 | A connectogram representing bilateral connectivity patterns of the ventral ILF subfascicle. Atlas regions are listed around the circumference of the connectogram. TInf, inferior temporal gyrus; TMd, middle temporal gyrus; TSp, superior temporal gyrus; Fu, fusiform gyrus; OInf, inferior occipital gyrus; OMd, middle occipital gyrus; OSp, superior occipital gyrus; Cal, calcarine gyrus; Cun, cuneus; Lin, lingual gyrus. Left and Right hemispheric connections are demonstrated with a suffix \_L or \_R.



Darker color represents stronger index. (NB, Left and Right hemispheric connections are demonstrated with a suffix \_L or \_R). TInf, inferior temporal gyrus; TMd, middle temporal gyrus; TSp, superior temporal gyrus; Fu, fusiform gyrus; OInf, inferior occipital gyrus; OMd, middle occipital gyrus; OSp, superior occipital gyrus; Cal, calcarine gyrus; Cun, cuneus; Lin, lingual gyrus; Am – amygdala; Hipp, hippocampus.

vertical part of the SLF, also known as the parietal aslant tract, was exposed and removed to further expose the temporal AF, which was subsequently removed to expose underlying ILF fibers. ILF fibers were exposed along their antero-posterior course to the occipital pole. Toward the occipital pole, ascending fibers belonging to the VOF, were removed. From an inferior aspect U-fibers from the fusiform gyrus were gently removed to expose longitudinal ILF fibers. Lateral and medial occipito-temporal sulci were opened, as well as the lingual and parahippocampal gyri. The ILF represented the most lateral and ventral tract among other fascicles within the sagittal stratum. Its bifurcation in the sagittal plane led to identification of two components: ventral and dorsal. Dissection results reinforced the tractographic results: The ILF lays ventro-laterally to the posterior sagittal stratum fibers, all coursing in parallel to terminate within the occipital lobe. The ILF ran lateral to the temporal horn and its floor. At the temporal pole, the ILF ended lateral to the UF. At the occipital pole, it was flattened between the VOF laterally, and the IFOF medially (**Figures 7C,D**).

### DISCUSSION

### ILF Morphology, Subdivisions and Relations

We have demonstrated that GQI-tractography can reproduce bilateral ILF bundles according to previous descriptions. Our results confirm the existence of a robust, temporo-occipital fascicle originating from medial and lateral anterior temporal areas and terminating throughout the occipital lobe. We did

not find any evidence of the 'indirect' U-fiber component of the ILF (Tusa and Ungerleider, 1985; Catani et al., 2003). Further, we were able to divide the whole bundles into smaller sub-fascicles for individual study. We used the clear posterior bifurcation of the ILF into ventral and dorsal terminations in the temporo-occipital cortex. Though the ventral sub-fascicle itself demonstrated a bi-pronged 'lateral-medial' termination pattern, we elected to simplify our approach using dorsal and ventral divisions for further analysis. For both left and right hemispheric ILFs, posterior arrangement was consistent over 30 subjects and within the HCP 1021 atlas. Anterior ILF termination pattern was relatively more variable and did not demonstrate consistent inter-subject morphology, permitting subdivision. In terms of morphology and subdivision, our results are generally consistent with both older (Catani et al., 2002, 2003; Fernández-Miranda et al., 2008) and more recent (Latini, 2015; Latini et al., 2017) studies, all of which demonstrate posterior terminations within both dorsal (lateral and medial occipital cortex) and ventral (i.e., lateral, temporo-occipital areas and medial, basal-occipital areas). Previous authors have elected to subdivide the ILF based upon connectivity profiles. For example, Catani et al. (2003), Latini (2015), and Latini et al. (2017) defined a 'cuneate' branch of the dorsal ILF. From our connectivity analysis, CI's were relatively small for cuneate connections of the left and right dorsal, ventral

anterior temporal lobe consists of fibers which are part of the ILF.

and undivided ILFs, compared to other connections. Using connective rather than morphological subdivisions introduces complexity into a nomenclature system, as an arbitrary threshold must be chosen to define what CI qualifies a connection as a significant 'subdivision' of the ILF. Regarding relations, we confirm that the ILF is the most superficial and ventral white matter fasciculus of the temporal lobe. Though cortico-cortical U-fibers may lie superficial to the robust ILF and indirectly connect temporal with occipital areas via short connections, we have found no tractographic evidence to suggest their anatomical affiliation with the ILF system.

### Quantitative ILF Connectivity

We used a purely quantitative method of assessing the connectivity of the ILF and its sub-fasciculi. This represents a methodological evolution from our previous studies (Wang et al., 2013, 2016; Fernández-Miranda et al., 2015; Panesar et al., 2017), all of which used qualitative correlations of fiber end-points within atlas regions. Not only does our quantitative technique give an objective measure of cortical termination, it also allows a metric analysis of the strength of connectivity between cortical areas. The need for quantitative termination analysis methods in tractographic studies has been identified (Girard et al., 2014) and will ultimately permit mapping of the connectome (Behrens

and Sporns, 2012). Recently published quantitative analysis methods include statistical analysis of tract termination patterns derived from whole brain probabilistic DTI tractograms (Hau et al., 2016, 2017). Our use of connectograms with deterministic tractography (Yeh et al., 2018) aids understanding of connectivity strength by visually demonstrating bilateral ILF temporooccipital connectivity in a topographically correlated matrix.

From our connectivity analysis of individual sub-fasciculi, it was clear that both sub-fascicles, and the undivided ILF were leftward-lateralized in terms of connectivity. This is confirmed by an overall increased average connectivity strength on the left versus the right for individual sub-fascicles. Though the differences in connectivity between left and right ILF subcomponents were not statistically significant, we postulate that this arose due to a relatively small size of the subject set and large spread of CIs. Moreover, though there were some fibers traversing between the cuneus, lingual gyri, and the amygdala, these connections were negligible in relative CI (see also section "ILF Morphology, Subdivisions and Relations"). We therefore question whether the 'Li-Am' bundle as proposed by Latini (2015) or 'cuneate' bundle as proposed by Catani et al. (2003) should be attributed as unique subdivisions of the ILF, versus simple connections.

Our connectivity results demonstrate the novel finding of strong leftward connectivity of the ventral ILF between the calcarine area and all three temporal gyri. This is in contrast to existing views (Catani et al., 2003; Latini, 2015; Latini et al., 2017). Further, connections originating from the calcarine area had strong relative CIs for the ventral sub-fascicle, and also when the undivided ILF was analyzed as a whole. The functional implications of these connections are discussed in Section "ILF Volumetry." Another novel finding is minimal or absent connectivity with the antero-medial temporal areas, which is again in contrast to previous postulations (Catani et al., 2003; Mills et al., 2013; Latini et al., 2017; Jennings et al., 2018). We postulate that this is because medial anterior-temporal white matter trajectories are occupied primarily by the uncinate fasciculus (**Figure 6D**). The ventral uncinate fasciculus descends via the temporal stem and occupies medial antero-temporal white matter trajectories (Panesar et al., 2017). As the ILF is the most

superficial, ventral white matter tract, its connectivity is primarily to the lateral cortical aspects of the temporal lobe. As all of these studies used DTI, we attribute previous findings of ILF anteromedial temporal connectivity to inability of DTI to accurately discern close proximity fiber systems at a high resolution versus our GQI method.

#### ILF Volumetry

Our volumetric results, demonstrating significant leftward asymmetry of the whole ILF, and its individual constituents are in concordance with earlier DTI-derived results (Wakana et al., 2007; Thiebaut de Schotten et al., 2011), however, it was in contrast to recent postulations by Latini et al. (2017), which demonstrated rightward volumetric asymmetry. Our results of leftward volumetric asymmetry can be attributed to overall increased leftward volume of both ILF subcomponents, which all demonstrate LIs >0. Our results are in concordance with our previous study into IFOF structure, which also demonstrated leftward lateralization of this ventral tract. The ILF and IFOF may therefore be both structurally and functionally related.

#### Functional Postulations

Regarding function of individual subfascicles, the left dorsal ILF demonstrated strongest connectivity between the superior occipital gyrus and the two dorsal temporal gyri. The superior temporal gyrus has been implicated in emotional responses to visual facial stimuli (Domínguez-Borràs et al., 2009). The middle temporal gyrus has been implicated in both semantic access (Davey et al., 2015), distance judgment (Vandenberghe et al., 1996), and facial recognition (Haxby et al., 2000). The dorsal ILF is therefore strong candidate for a neural substrate subserving both semantic, spatial and facial recognition tasks. The left ventral sub-fascicle connected strongly with inferior and middle occipital gyri, cuneate, calcarine, and lingual areas, respectively. Both had extensive terminations within the three lateral temporal gyri. Right sided ventral ILFs demonstrated exclusive connectivity between lingual and the all three temporal gyri. The cuneus is implicated in extrastriate visual processing (Allison et al., 1994). The calcarine cortex contains V1, and the lingual gyrus has been implicated in both facial recognition (Haxby et al., 2000) on the right and binocular spatial orientation bilaterally (Lee et al., 2001). Our results therefore implicate bilateral ventral components of the ILF to be implicated in facial recognition and would explain lesion-induced prosopagnosia (Benson, 1974; Bauer, 1984; Michel et al., 1989; Takahashi et al., 1995; Grossi et al., 2014).

Within the 'dorsal-ventral' model of language functionality, prevailing opinion is of a leftward-lateralized ventral semantic system (Hickok and Poeppel, 2015). Our connective and volumetric results therefore support a proposal of the ILFs role within the ventral semantic-stream. Further, due to its rich occipito-temporal connectivity, our results also support a dominant ILF role as an integrative system between the visual cortices and the relevant lateral temporal gyri, as evidenced by lesional studies (Bonilha et al., 2017). Recent theories (Dehaene and Cohen, 2011; Zemmoura et al., 2015) have proposed the leftsided occipital visual word form area as a 'functionally recycled' neural substrate, originally evolved to subserve facial recognition. This theory is consistent with theories of ILF function, which implicate it in both semantic language and visual recognition roles. Our findings of ILF leftward-lateralization also reinforce this postulation.

### Technical Considerations and Limitations

Our method was intended to generate whole ILFs without a priori regarding either connectivity or subdivision. Subsequent separation of the bundles did indeed reveal a dorsal-ventral ILF arrangement in accordance with the literature. Nevertheless, based upon our connectivity results, nomenclature for subdivided sub-fasciculi should not be used to indicate connectivity patterns (see sections "ILF Morphology, Subdivisions and Relations," "Quantitative ILF Connectivity"), but instead only morphological divisions of the ILF. This is exemplified by the right ventral ILF sub-fascicles, both of which demonstrate connectivity to both medial and lateral occipital cortices. Regarding connectivity overlaps between the divisions, whether these occur due to true segregation between sub-fasciculi or whether these are due to tractographic artifacts remains to be determined. In an effort to reconcile this issue, and also proposed by Latini et al. (2017), we emphasize that both connectivity and volumetric results for whole, merged ILFs take precedent over results for individual sub-tract analyses. Though connective segmentation of the ILF may seem appealing, especially using our CI metric, we opt not to use it to define specific 'subdivisions' of the ILF, as we cannot define the arbitrary threshold deeming a connection significant enough to be considered an independent sub-fasciculus rather than a common connection or artifact.

Our current high-resolution GQI based imaging technique gives ability to accurately track close-proximity and crossing fiber systems to termination points within cortical mantle. The combination of this technique with quantitative connectivity analyses is the basis for connectogram construction and provides our most detailed and objective study to date. Further, when the pitfalls of DTI imaging are taken into consideration, the combination of tractography with objective connectivity analyses offers a superior solution for in vivo research studies. Post-mortem fiber dissection has played a pivotal role in neuroanatomical studies and remains a valuable technique for comparison and validation of tractographic findings. Nonetheless, post-mortem fiber dissection provides only qualitative evaluation and is prone to human errors. It is becoming rapidly surpassed by tractographic methods that allow for quantifying neuroanatomic parameters such as connectivity, and within large subject-sets as we have demonstrated.

#### CONCLUSION

We have successfully replicated a subdivided ILF bundle consistently over the range of 30 subjects. Though our results bear

resemblance to previous postulations regarding the subdivision and gross connectivity of the ILF, we have offered a detailed picture of its bilateral connectivity patterns. The ILF is a connectively and volumetrically leftward-lateralized, ventral temporal bundle with rich occipito-temporal connectivity. Though tractography cannot give direct functional insight, our connectivity analysis taken in context with functional data pertaining to cortical functions offer an accurate and confirmatory insight into its postulated roles. We conclude by calling for increased adoption of GQI tractography and quantitative connectometry for anatomical studies.

#### ETHICS STATEMENT

This study was carried out in accordance with the recommendations and approval of the Institutional Review Board

#### REFERENCES


at the University of Pittsburgh. All subjects gave written informed consent in accordance with the Declaration of Helsinki.

#### AUTHOR CONTRIBUTIONS

SP contributed to the main data collection, writing, image creation, and tabulation. F-CY designed the algorithm, calculated the connectivity indices, and created the connectogram. TJ contributed to the WM dissection and pictures and discussion. WH contributed to the discussion. JF-M was the PI of the lab and contributed the main ideas, to document editing, and finalization.

### FUNDING

This study was funded by the Pittsburgh Brain Institute.

to move beyond DTI. J. Neurosurg. 118, 1367–1377. doi: 10.3171/2013.2. JNS121294


physics, fiber tract anatomy, and tumor imaging patterns. Am. J. Neuroradiol. 25, 356–369.


evidence from diffusion tensor imaging. Hum. Brain Mapp. 34, 1025–1034. doi: 10.1002/hbm.21492


**Conflict of Interest Statement:** 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.

Copyright © 2018 Panesar, Yeh, Jacquesson, Hula and Fernandez-Miranda. 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 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.

# Derivation of Fiber Orientations From Oblique Views Through Human Brain Sections in 3D-Polarized Light Imaging

Daniel Schmitz <sup>1</sup> \*, Sascha E. A. Muenzing<sup>1</sup> , Martin Schober <sup>1</sup> , Nicole Schubert <sup>1</sup> , Martina Minnerop1,2, Thomas Lippert 3,4, Katrin Amunts 1,5 and Markus Axer <sup>1</sup>

1 Institute of Neuroscience and Medicine-1 (INM-1), Forschungszentrum Jülich, Jülich, Germany, <sup>2</sup> Center for Movement Disorders and Neuromodulation, Department of Neurology and Institute of Clinical Neuroscience and Medical Psychology, Medical Faculty, Heinrich-Heine University, Düsseldorf, Germany, <sup>3</sup> Jülich Supercomputing Center, Forschungszentrum Jülich, Jülich, Germany, <sup>4</sup> Bergische Universität Wuppertal, Wuppertal, Germany, <sup>5</sup> C. and O. Vogt Institute for Brain Research, Medical Faculty, Heinrich-Heine University Düsseldorf, Düsseldorf, Germany

3D-Polarized Light Imaging (3D-PLI) enables high-resolution three-dimensional mapping of the nerve fiber architecture in unstained histological brain sections based on the intrinsic birefringence of myelinated nerve fibers. The interpretation of the measured birefringent signals comes with conjointly measured information about the local fiber birefringence strength and the fiber orientation. In this study, we present a novel approach to disentangle both parameters from each other based on a weighted least squares routine (ROFL) applied to oblique polarimetric 3D-PLI measurements. This approach was compared to a previously described analytical method on simulated and experimental data obtained from a post mortem human brain. Analysis of the simulations revealed in case of ROFL a distinctly increased level of confidence to determine steep and flat fiber orientations with respect to the brain sectioning plane. Based on analysis of histological sections of a human brain dataset, it was demonstrated that ROFL provides a coherent characterization of cortical, subcortical, and white matter regions in terms of fiber orientation and birefringence strength, within and across sections. Oblique measurements combined with ROFL analysis opens up new ways to determine physical brain tissue properties by means of 3D-PLI microscopy.

Keywords: neuroimaging, modeling, 3D-PLI, white matter anatomy, fiber architecture

## 1. INTRODUCTION

Understanding the human brain's function and dysfunction requires a thorough knowledge about the brain's fiber tracts, forming a dense network of connections within, but also between the different brain regions. Over the last years, several imaging techniques have emerged which are capable of resolving anatomical structures with different spatial resolutions. At the millimeter scale, diffusion MRI is the most prominent one as it is applicable to both in vivo and post mortem brains (Basser et al., 1994; Johansen-Berg and Behrens, 2009; Mori and Tournier, 2014). At ultra-high resolution two-photon microscopy (Laperchia et al., 2013), light-sheet microscopy (Silvestri et al., 2012), and electron microscopy (Knott et al., 2008), amongst others, have been exploited to image single cells and neurons in 3D space as well as their local connections. Yet these techniques come with the cost of excessive measurement time, large amounts of data and limited fields of view (lateral and axial), impeding the study of larger brain volumes so far.

#### Edited by:

Laurent Petit, Centre National de la Recherche Scientifique (CNRS), France

#### Reviewed by:

Karla Miller, University of Oxford, United Kingdom Caroline Magnain, Massachusetts General Hospital, Harvard Medical School, United States

> \*Correspondence: Daniel Schmitz da.schmitz@fz-juelich.de

Received: 07 May 2018 Accepted: 27 August 2018 Published: 27 September 2018

#### Citation:

Schmitz D, Muenzing SEA, Schober M, Schubert N, Minnerop M, Lippert T, Amunts K and Axer M (2018) Derivation of Fiber Orientations From Oblique Views Through Human Brain Sections in 3D-Polarized Light Imaging. Front. Neuroanat. 12:75. doi: 10.3389/fnana.2018.00075

3D-Polarized Light Imaging (3D-PLI) (Axer et al., 2011a,b) is a microscopic technique that bridges the gap between largescale imaging techniques such as diffusion MRI and ultra-high resolution microscopy. It enables the reconstruction of fiber tracts at the meso- and micro-scale from unstained histologial sections. Recently, polarization sensitive optical coherence tomography (PSOCT) has emerged as a promising technique for the mapping of fiber bundles with the ability of depth-resolved scanning of brain blocks (Wang et al., 2011, 2014; Magnain et al., 2014). While the signals measured by PSOCT and 3D-PLI both arise due to the birefringence of brain tissue, their fundamental difference is that PSOCT captures the reflected light instead of the transmitted light as 3D-PLI (Caspers and Axer, 2017).

A pitfall of any histological imaging technique is the extraction of information in the direction perpendicular to the sectioning plane, i.e., in the depths of the histological section. While inplane fiber orientations can be obtained from texture information of microscopic images via the structure tensor without the need for any biophysical model (Budde and Frank, 2012), little information is directly available about orientations perpendicular to the sectioning plane. Lately, three-dimensional structure analysis has been applied to image stacks obtained from light microscopy (Khan et al., 2015) and confocal microscopy (Schilling et al., 2016). This comes with several difficulties: the prerequisites of a precise 3D-alignment of the image stack as well as a coherent contrast over the whole image stack and the strong dependency of the obtained orientations on the chosen pre-processing pipeline and kernel size of the structure tensor (Khan et al., 2015).

As an alternative, 3D-PLI enables to determine fiber orientations directly from the measurement data without the need for a manual choice of image processing parameters. This is achieved by taking measurement data from different views by tilting the brain section and analyzing it with a biophysical model (Axer et al., 2011a). Especially, the estimation of the outof-sectioning-plane orientation (inclination angle) benefits from the additional information gained by tilting the brain section as it is measured conjointly with the fiber density.

Several algorithms for the analysis of the data acquired with the tiltable specimen stage have been presented so far. The Markov Random Field algorithm by Kleiner et al. (2012) and the total variation-cased method by Alimi et al. (2017) both resulted in a robust estimation of the inclination sign. While this solved the inclination sign ambiguity, the absolute inclination was still undetermined. This problem was first solved by Wiese et al. (2014) who showed how the tilted measurements can be utilized to extract inclination and fiber density independently from each other. Their algorithm is based on a discrete Fourier transformation of the different tilting positions (denoted as DFT algorithm). Due to its analytical nature, the DFT algorithm is computationally efficient but suffers from noise instability.

The here presented approach seeks to overcome this noise instability and provide a more robust separation of fiber density and orientation. Our approach utilizes a weighted least squares algorithm to process the tilted measurements. It is then compared to the analytical DFT algorithm on synthetic and experimental data. The examined experimental results represent the first analysis of large-scale 3D-reconstructed human brain data in 3D-PLI. The new approach was introduced in Schmitz et al. (2018) for the first time, yet this work represents a vast extension of the analysis and results.

## 2. MATERIALS AND METHODS

## 2.1. Least Squares Estimation of Fiber Parameters

#### 2.1.1. 3D-PLI

3D-PLI utilizes the birefringence of nerve fibers which is measured in customized polarimeters. The birefringence originates from the regular arrangement of lipids in the myelin sheath resulting in optical anisotropy. This anisotropy causes a phase shift of incident polarized light passing the brain tissue. The optical setup as described in Axer et al. (2011b) is depicted in **Figure 1**: first unpolarized light from an LED array (custommade design, FZJ-SSQ300-ALK-G, iiM, Germany) passes a first polarization filter and a quarter-wave retarder mounted with a principle axis offset of 45◦ with respect to the polarizer (Jos. Schneider Optische Werke GmbH, Germany), circularly polarizing the light. Then the light interacts with the brain tissue mounted on a tilting stage. A second polarizer, rotated 90◦ with respect to the first one serves as analyzer. The outgoing light is captured with a CCD camera with 14 bit depth (AxioCam HRc, Zeiss, Germany), capturing images at a pixel size of 64 × 64µm. The filters and the retarder are rotated simultaneously in steps of ρ = 10◦ , yielding an image series acquired by the camera.

The effective physical model behind 3D-PLI describes myelinated nerve fibers as negative uniaxial birefringent crystals. Assuming ideal optical components, this model yields the following expression of the light intensity in each pixel during rotation of the filters (Axer et al., 2011a):

$$I(\rho,\varphi,\delta) = \frac{I\_T}{2} \cdot \left(1 + \sin(2(\rho - \varphi)) \cdot \sin\delta\right) \tag{1}$$

with the transmittance <sup>I</sup><sup>T</sup> 2 , the in-plane orientation ϕ (direction angle), and the retardation δ. For the retardation the following expression is presumed:

$$\delta = \frac{\pi}{2} d \cos^2(\alpha) \tag{2}$$

d with the relative thickness d and the out-of-plane orientation α (inclination angle). For readability purposes we denote d as trel in contrast to Axer et al. (2011a). The angles are defined in the range ϕ ∈ [0◦ , 180◦ ] and α ∈ [−90◦ , 90◦ ]. d was introduced by Axer et al. (2011a) as the combined effect of section thickness t<sup>s</sup> , birefringence 1n, and illumination wavelength λ:

$$d = 4 \frac{\mathfrak{t}\_s \Delta n}{\lambda} \tag{3}$$

The sinusoidal profile is analyzed with a Fourier analysis: the I<sup>T</sup> is given by the average of the profile, the direction angle by its phase

and the retardation sin δ by its relative amplitude (Axer et al., 2011a). While this analysis enables unambiguous determination of transmittance, direction angle, and retardation, α and d are mapped onto one value, the retardation. Disentanglement of both parameters is a priori not possible, as any combination of d and α which fulfills Equation (2) is valid solution. Furthermore, for d > 1 the outer sin function induces an additional ambiguity. Therefore, the separation of d and α is only bijective for d ∈ (0, 1]. In Axer et al. (2011a) a constant relative thickness which was extracted by statistical approach is assumed. This enables the calculation of the inclination by inverting Equation (2). While this approach offers a good first guess for the inclinations, the assumption of a constant d over the whole section did not take local differences between distinct brain regions into account. For example, the cortex contains far less myelinated axons than dense white matter regions such as the corpus callosum. Also, it has to be noted that the relative thickness depends on the section thickness which is not precisely known and the birefringence which depends on the local amount of myelin (Goethlin, 1913). Therefore, it is clear that the unbiased disentanglement of both parameters requires additional measurement information.

A mechanical solution to gain more measurement information is a tiltable specimen stage built into the polarimetric system, which introduces a pre-defined tilting angle to the nerve fiber orientation relative to the sectioning plane (Axer et al., 2011a). By this means, polarimetric measurements can be performed from oblique views. Experimentally, this is realized by a two-nested axis system enabling rotations of the specimen around the x- and y-axis by up to 8◦ . In routine 4 tilting positions are measured: at each position the brain section is tilted by ±8 ◦ with respect to one of the axes, denoted as N(orth), S(outh), E(ast), and W(est) (see **Figure 1**). All measurements taken with a tilted sample are registered on the measurement without tilt (planar measurement) by a projective linear transformation.

#### 2.1.2. Light Intensity Distribution of the Imaging System

For an accurate model of the data acquisition, the noise level of the imaging system must be taken into account. As the detected signals are light intensities captured by a CCD camera in our case, the dominant noise sources are shot noise during image acquisition and the internal signal processing of the camera (Bertolotti, 1974; Goodman, 2000). Here, the distribution of interest is the distribution of the number of detected photoelectrons during exposure time per pixel. With increasing photon count the variance of the photon count increases (Goodman, 2000). The dependency of the variance of the detected photons σ <sup>2</sup> on the number of detected photons as the expected value µ can only be determined experimentally. Therefore, 100 samples of the same scene were taken and analyzed pixelwise for the variances and expected values of the measured light intensities. The resulting relationship between variance and expected value is given by σ <sup>2</sup> = 3µ (cf. **Supplementary Material** for details). The multiplicative factor which relates variance and expected value is called gain factor g (σ <sup>2</sup> = gµ), in our case g = 3. A common way to model overdispersed count data is a negative binomial model as it allows an arbitrary expectation value and variance. In general, the occurrence of k under a negative binomial distribution parameterized by the expected value and variance is given by

$$P(k|\mu,\sigma) = \binom{k-1+\frac{\mu^2}{\sigma^2-\mu}}{k} \left(\frac{\sigma^2-\mu}{\sigma^2}\right)^k \left(\frac{\mu}{\sigma^2}\right)^{\frac{\mu^2}{\sigma^2-\mu}} \tag{4}$$

The probability to observe the light intensity I given the expected light intensity µ and the gain factor g can then be expressed as

$$P(I|\mu,\mathcal{g}) = \binom{I-1+\frac{\mu}{\mathcal{g}-1}}{I} \left(\frac{\mathcal{g}-1}{\mathcal{g}}\right)^I \mathcal{g}^{-\frac{\mu}{\mathcal{g}-1}}\tag{5}$$

#### 2.1.3. The Tilting Coordinate System

We will use the coordinate system introduced by Kleiner et al. (2012): it describes the tilting stage by the direction in which the stage is tilted ψ and the actual tilting angle τ as sketched in **Figure 1**. For the four tilting directions carried out in the standard measurement, the tilting directions then are ψ = 0, 90, 180, 270◦ with a tilting angle of τ = 8 ◦ . With NTilt as the number of tilting positions the tilting directions equidistantly spanning the space are in general given by ψ<sup>j</sup> = 2π NTilt (j − 1), j ∈ [1, 2, ... , NTilt] with the index j indicating the tilting position. The tilted vector **r<sup>t</sup>** is obtained by applying the appropriate rotation to the vector in the planar position **r**: **r<sup>t</sup>** = **R**(ψ, τ )**r**. The full rotation matrix **R**(ψ, τ ) is derived by first rotating around the z-axis by −ψ, then rotating around the y-axis by τ and then rotating back around the z-axis by ψ (Wiese, 2017):

**R**(ψ, τ ) = **R** z (ψ)**R** y (τ )**R** z (−ψ)

$$\begin{aligned} \mathbf{x} &= \begin{pmatrix} \cos(\psi) & -\sin(\psi) & 0 \\ \sin(\psi) & \cos(\psi) & 0 \\ 0 & 0 & 1 \end{pmatrix} \cdot \begin{pmatrix} \cos(\mathbf{r}) & 0 & \sin(\mathbf{r}) \\ 0 & 1 & 0 \\ -\sin(\mathbf{r}) & 0 & \cos(\mathbf{r}) \end{pmatrix} \cdot \begin{pmatrix} \cos(\psi) & \sin(\psi) & 0 \\ -\sin(\psi) & \cos(\psi) & 0 \\ 0 & 0 & 1 \end{pmatrix} \\ &= \begin{pmatrix} \cos(\mathbf{r})\cos(\psi)^2 + \sin(\psi)^2 & (\cos(\mathbf{r}) - 1)\sin(\psi)\cos(\psi) & \cos(\psi)\sin(\mathbf{r}) \\ (\cos(\mathbf{r}) - 1)\sin(\psi)\cos(\psi) & \cos(\mathbf{r})\sin(\psi)^2 + \cos(\psi)^2 & \sin(\psi)\sin(\mathbf{r}) \\ -\cos(\psi)\sin(\mathbf{r}) & -\sin(\psi)\sin(\mathbf{r}) & \cos(\mathbf{r}) \end{pmatrix} . \end{aligned} \tag{6}$$

As stated in Wiese et al. (2014), due to refraction at the brain tissue, the actual tilting angle in the sample is reduced and has to be adjusted according to Snell's law. Assuming a refractive index of n ≈ 1.45 for human tissue based on studies of de Campos Vidal et al. (1980), the tilting angle of the light in the tissue is given by <sup>τ</sup>int <sup>≈</sup> 5.51◦ for a tilt of the tilting stage by τ = 8 ◦ .

#### 2.1.4. Least-Squares Algorithm

We introduce the indices j for the tilting direction (including 0 for the planar measurement) and i for the rotation angle of the polarization filters. Additionaly, we denote the total number of measured tilting stage positions (tilted measurements and the planar measurement) as N<sup>T</sup> = NTilt + 1 and the number of polarizer positions acquired as NP. All measurement data accumulates to N<sup>T</sup> · N<sup>P</sup> light intensities Iji for each pixel. Application of the Fourier analysis to all measurements as described in Axer et al. (2011a) results in N<sup>T</sup> transmittance values Ij,<sup>T</sup> , N<sup>T</sup> retardation values |sin δ<sup>j</sup> | and N<sup>T</sup> direction values ϕ<sup>j</sup> . In this notation the intensity curve of a measurement can than be expressed as

$$I\_{ji}(\rho\_i, \varphi\_j, \alpha\_j, d\_j) = \frac{I\_{j,T}}{2} \left( 1 + \sin(2(\rho\_i - \varphi\_j)) \sin\left(\frac{\pi}{2} d\_j \cos^2(\alpha\_j)\right) \right) \tag{7}$$

In a first step, the effect of the average light intensity, the transmittance, is eliminated from the fiber orientation determination. Eliminating it from the further process is necessary as the light experiences additional absorption and refraction effects in a tilted measurement, which cannot be included in the model. Therefore, we define the normalized light intensity

$$I\_{N\_{ji}} = \frac{2I\_{ji}}{I\_{j,T}} - 1 \tag{8}$$

which shifts the data into the range [−1, 1]. As the variance of the light intensities is known as σ 2 Iji = gIji the variance of the normalized light intensity can be derived using error propagation as

$$
\sigma\_{N\_{ji}}^2 = \frac{gI\_{ji}}{I\_{j,T}^2} + \frac{gI\_{ji}^2}{N\_P I\_{j,T}^3} \tag{9}
$$

The normalized light intensities fji predicted by the 3D-PLI model are given by inserting Equation (7) into Equation (8):

$$f\_{ji}(\varphi\_j, \alpha\_j, d\_j, \rho\_i) = \sin(2 \cdot (\rho\_i - \varphi\_j)) \sin\left(\frac{\pi}{2} \cdot d\_j \cdot \cos(\alpha\_j)^2\right) \tag{10}$$

Fitting f to the normalized light intensities I<sup>N</sup> can now be formulated as a weighted least squares problem. With the weights wij = σ −1 Nji the optimization problem is given by

$$\underset{\varphi,\alpha,d}{\text{argmin}}\,\chi^2 = \underset{\varphi,\alpha,d}{\text{argmin}}\sum\_{j=0}^{N\_T}\sum\_{i=0}^{N\_P} \left( \left( f\_{ji}(\varphi\_j,\alpha\_j,d\_j,\rho\_i) - I\_{N\_{ji}} \right) \cdot \boldsymbol{\nu}\_{ji} \right)^2 \tag{11}$$

subject to ϕ ∈ [0, π], α ∈ [− π 2 , π 2 ], d ≥ 0. At this point we do not restrict the relative thickness to d ≤ 1 as a value of d > 1 could also occur if the current 3D-PLI model simply does not describe the data well enough. Two issues need to be overcome to solve the optimization problem: the bounded parameter space and a suitable starting point for the local optimization. These issues resemble the optimization problem of the maximum likelihood algorithm described in Wiese (2017).

A starting point for the direction angle is given by the direction angle derived from the planar measurement ϕ0. For the inclination α and the relative thickness d, the starting point is determined by brute force minimization. Based on simulation studies, a 6 × 6 grid equidistantly spanning the parameter space ([ϕ0, α<sup>l</sup> , d<sup>k</sup> ], k, l = 1, ... , 6) was found to result in a promising first guess for the subsequent optimization.

The boundaries of the parameter space could be accounted for by utilizing an optimization algorithm capable of dealing with hard boundaries. Yet, in our case a more elegant solution is to exploit the symmetry of the problem. This enables to reformulate the bounded optimization problem as an unbounded problem. As f(ϕ, α, d) = −f(ϕ, α, −d), it is sufficient to allow all relative thicknesses d and take the absolute value of the relative section thickness d = |d| before calculating the normalized intensities. For the fiber orientations, it is also possible to allow all orientations considering the symmetry of spherical coordinates. The unbounded orientation given as (ϕu, αu) can be transformed back into the standard 3D-PLI parameter space by the following transformations:

$$\alpha = \left( \left( \alpha\_{\iota} + \frac{\pi}{2} \right) \bmod \pi - \frac{\pi}{2} \right) \text{sgn} \left( \frac{1}{2} - \left\lfloor \frac{\varphi\_{\iota}}{\pi} \bmod 2 \right) \right) \tag{12}$$

$$
\varphi = \varphi\_{\mu} \bmod \pi \tag{13}
$$

Before calculating the normalized light intensities f , these transformations are used to symmetrize the unbounded **Algorithm 1:** Pseudocode of the ROFL algorithm

**Data**: Maps of ϕ0, Ij,T, Iji , j ∈ {0, ... , NT}, i ∈ {0, ... , NP} **Result**: Maps of ϕ, α, d, χ 2 **for** every image pixel **do** // calculate normalized intensities and their standard deviations **for** j = 0 : N<sup>T</sup> **do for** i = 0 : N<sup>P</sup> **do** INji = Iji 2Ij,<sup>T</sup> − 1 σNji = error prop. of σIj,<sup>i</sup> , σIj,<sup>T</sup> to INj,<sup>i</sup> // brute force minimization **for** α<sup>l</sup> , d<sup>k</sup> ∈ bruteforce grid **do** χ 2 <sup>l</sup>,<sup>k</sup> <sup>=</sup> <sup>χ</sup> 2 ϕ0, α<sup>l</sup> , d<sup>k</sup> , (IN<sup>00</sup> , σN<sup>00</sup> , . ..) α0, d<sup>0</sup> = argmin α,d χ 2 l,k // optimization Levenberg-Marquardt optimization of χ <sup>2</sup> with initial point (ϕ0, α0, d0) → ϕ, α, d, χ 2 transform ϕ, α, d into 3D-PLI parameter space

orientation into the parameter space. As optimization algorithm the Levenberg-Marquardt algorithm was chosen (Levenberg, 1944; Marquardt, 1963). In case of convergence to a point outside of the parameter space, the before mentioned transformations are used to symmetrize the result back into the parameter space. This newly developped algorithm is denoted as Robust Orientation Fitting via Least Squares (ROFL). A pseudocode is given in **Algorithm 1**.

For this study, the ROFL algorithm was implemented in Python based on the packages NumPy (van der Walt et al., 2011) and SciPy (Jones et al., 2001 ). The necessary computation time is very high as the optimization has to be carried out in every image pixel (about 4 core hours for an image size of ≈ 10<sup>6</sup> pixels). Therefore, the algorithm was parallelized pixelwise using mpi4py (Dalcin et al., 2005) and all calculations were conducted on the JURECA system (Jülich Supercomputing Centre, 2016).

#### 2.2. Simulated Data

Monte-Carlo simulations were carried out to test the robustness of the presented approach and the DFT algorithm against measurement noise.

#### 2.2.1. Generation of Synthetic Data

Simulations can provide information about two problems. Firstly, the accuracy of the obtained parameters can be compared to the synthetic ground truth. Secondly, the algorithms can be tested against biases in their parameter estimation. One sythetic dataset each was created to answer these questions. For both datasets, the generation of synthetic 3D-PLI signals for a specific parameter set is executed in the same way.

Synthetic signals were generated in a similar manner as in Wiese et al. (2014). A ground truth orientation vector **v**<sup>g</sup> with the parameters (d<sup>g</sup> , ϕ<sup>g</sup> , α<sup>g</sup> ) was rotated into four tilting positions (North, South, East, West) by τ = 5.51◦ , the assumed internal tilting angle for human brain tissue. For each tilting position a sinusoidal light intensity profile according to Equation (1) was calculated. For the transmittance I<sup>T</sup> 2500 was chosen, as it represents a typical transmittance value for human brain sections. To mimic the experimentally observed light intensity distribution, for each calculated light intensity I a noisy light intensity Inoisy was then computed by drawing one sample from a negative binomial distribution with expected value given by I and variance 3I: Inoisy ∝ N B(I, 3I).

In a first simulation, the accuracy of the obtained parameters was assessed. As our primary interest here was the validation of the reconstruction of the inclination and the relative section thickness, only one direction angle, ϕ = 45◦ , was simulated. This direction angle also represents the worst case scenario for the four tilting positions as the angle between the orientation vector and the rotated orientation vector is maximal for ϕ = ψ and decreases with |ϕ − ψ|. With respect to the inclination angle, it is sufficient to simulate only positive inclinations as the inclination sign does not effect the reconstruction precision. As ground truth vectors all combinations of the fixed direction angle 45◦ , inclinations α from 0◦ , 1◦ , ... , 89◦ , 90◦ , and thicknesses d from 0.01, 0.02, ... , 0.89, 0.9 were simulated. These fiber configurations display all different inclination angles for different tissue scenarios characterized by the relative thickness d. For each ground truth vector 100,000 samples of 3D-PLI signals were generated utilizing the beforementioned method to enable a statistical analysis.

To tests the algorithms against biases in their inclination estimation, a second dataset of 500,000 samples of orientation vectors uniformly distributed on the unit sphere were computed. From these vectors, a dataset of direction and inclination angles was derived and used to calculate the noisy sinusoidal light intensity profiles with a relative thickness of d = 0.5.

#### 2.2.2. Evaluation

The simulated sinusoidal profiles were analyzed with the ROFL algorithm and the DFT algorithm resulting in the reconstructed vector **v**<sup>r</sup> with the parameters (d<sup>r</sup> , ϕ<sup>r</sup> , αr). The accuracy of the obtained orientation was then measured by the acute angle γ between the ground truth vector **v**<sup>g</sup> and the reconstructed vector **v**<sup>r</sup> :

$$\boldsymbol{\gamma} = \arccos(|\mathbf{v}\_r \cdot \mathbf{v}\_\xi|) \tag{14}$$

The acute angle respects the symmetry of the problem: as the parameter space is bounded to a half sphere, the maximal angle between two vectors is 90◦ . Therefore, the absolute value of the scalar product is used in Equation (14). For each ground truth vector, the overall reconstruction error is given by the mean deviation angle hγ i of all 100,000 samples. The reconstruction accuracy of the relative thickness d is evaluated by the absolute relative error between the obtained values and the ground truth: <sup>σ</sup><sup>d</sup> = | <sup>d</sup>r−d<sup>g</sup> dg |. The overall error for a fiber configuration is then again given by the mean relative error of all 100,000 samples hσdi.

To assess if the obtained inclinations are biased, inclination angles were reconstructed from the second synthetic dataset. The resulting inclination angle histograms were then compared to the ground truth inclination histogram.

### 2.3. Experimental Data

To test and compare the DFT and the ROFL algorithm on experimental data, a series of coronal sections of a right human hemisphere was investigated. The post mortem human tissue sample used for this study was acquired in accordance with the local ethic committee of our partner university at the Heinrich Heine University Düsseldorf. As confirmed by the ethic committee, postmortem human brain studies do not need any additional approval, if a written informed consent of the subject is available. For the research carried out here, this consent is available.

#### 2.3.1. Data Acquisition

The examined brain was removed within 24 h after the donor's death. The right hemisphere was fixed in 4% buffered formaldehyde solution for 15 days to prevent tissue degeneration. After immersion in a 20% solution of glycerin with Dimethyl sulfoxide (DMSO) for cryoprotection the brain was frozen at a temperature of −80◦C. The coronal cutting plane was decided to be orthogonal to a line connecting the anterior and posterior commissures. Before sectioning, the frontal lobe was cut off to create a stable cutting platform. The sectioning resulted in 843 sections of 70 µm thickness (Polycut CM 3500, Leica, Germany). Before cutting an image of the cutting surface (blockface image) was taken for each section. Each section was then mounted on a glass slide, immersed in a 20% solution of glycerin to avoid dehydration and sealed by a coverslip. All sections were measured with the mentioned setup with the standard measurement protocol of the planar measurement and four tilting positions.

#### 2.3.2. Data Analysis

After an intensity based calibration as described in Dammers et al. (2010) the Fourier analysis was utilized to obtain the standard 3D-PLI modalities transmittance, retardation and direction. Then the ROFL algorithm was used to extract maps of the direction and inclination angles, the relative section thickness d and the residuum χ 2 . As a comparison, the DFT algorithm was applied to determine the inclination angle and relative section thickness. Also, the residuum was calculated for the parameters obtained by DFT.

After the sectionwise analysis, the 3D-reconstruction of the measured sections requires a multi-step registration to retain the 3D volume of the measured brain. Here the blockface volume generated by aligning the blockface images to a 3D volume serves as a reference for the registration of the histological sections (Schober et al., 2015). In a first step, the histological sections are registered onto their corresponding blockface images via an affine transformation utilizing in-house developed software based on the software packages ITK (Yoo et al., 2002) and elastix (Klein et al., 2009). Then non-linear registration using the ANTs toolkit (Avants et al., 2011) was performed. Finally, the scalar modalities were spatially transformed using the obtained deformation fields. The direction angles were rotated according to the pixelwise rotations induced by the deformation computed by the nonlinear image registration. Pixel values were linearly interpolated. The reorientation can result in orientations lying outside of the parameter space, therefore the same mapping as in ROFL (see Equation 13) was applied after the reorientation.

#### 2.3.3. Validation

Traditional anatomical studies have described the fiber pathways globally and without considering differences between brains. At the level of detail presented in this study with a voxel size of 64 × 64 × 70µm<sup>3</sup> , the inter-subject variability in the fine structure of fiber orientations becomes much more relevant. Therefore, the resulting fiber orientations cannot easily be compared to anatomical atlases. A complementary measurement at the same resolution enabling a comparison is not available either. Additionally, a phantom for 3D-PLI has not been developed yet, impeding the possibility of a measurement with a known ground truth. Therefore, we chose to validate the resulting orientations based on their coherence across the whole volume and their alignment to anatomical boundaries, for example, within the sagittal stratum. The results of both algorithms were furthermore compared from a statistical point of view by the difference between the predicted and measured light intensities measured by the sum of the squared residuals.

## 3. RESULTS

### 3.1. Simulation

#### 3.1.1. Accuracy Evaluation

The simulation results of the first simulated dataset are depicted in **Figure 2**: The plots show the orientation reconstruction error and the relative reconstruction error of the relative thickness of both algorithms as a function of inclination angle α and relative thickness d.

#### **3.1.1.1. Orientation reconstruction**

The orientation reconstruction accuracy of the DFT is valley shaped with respect to the inclination angle and becomes minimal for an inclination angle of α ≈ 60◦ . From this minimum the orientation reconstruction error increases slightly for flat fibers. The most challenging orientations to analyze are very steep fibers with respect to the sectioning plane: for these, the orientations predicted by the DFT algorithm differ strongly from the ground truth. Even for relative thicknesses of d > 0.2 an average reconstruction error of hγ i ≈ 22◦ occurs for α > 80◦ . With decreasing relative thickness d the reconstruction accuracy also decreases: for d < 0.05, the obtained orientations are basically random.

In contrast, the orientation reconstruction error for the ROFL algorithm does not take the form of a valley: for all section thicknesses d > 0.06, the error is minimal for in-plane fibers with α = 0 ◦ and increases with the inclination angle. For very steep fibers of α > 80◦ , the error increases to γ ≈ 12◦ on average. As before, the reconstruction error increases with decreasing relative thickness d. In a direct comparison, the orientation reconstruction error achieved by the ROFL algorithm is lower

than the orientation reconstruction error achieved by the DFT algorithm for all configurations. The minimal error for the ROFL algorithm is given by hγ i = 1 ◦ , compared to hγ i = 1.6◦ for the DFT algorithm. For d ∈ [0.2, 0.9] and α ∈ [0◦ , 80◦ ], the orientation reconstruction accuracy achieved by the algorithms are on average 2◦ for ROFL and 4.5◦ for DFT with maximal values of 9.5◦ for ROFL and 18.8◦ for DFT.

#### **3.1.1.2. Reconstruction of the relative section thickness**

Similar to the the orientation reconstruction error, the relative reconstruction error hσdi of the relative thickness d increases with the inclination angle. Here, the errors obtained from the DFT algorithm do not assume the shape of a valley: besides a noise instability for very high relative thicknesses d > 0.8 the error increases with the inclination and increases with decreasing relative thickness d. The relative reconstruction error of the ROFL algorithm assumes the same shape but is lower: on average, the relative thickness of fibers with d ∈ [0.2, 0.9] and α ∈ [20◦ , 90◦ ] is determined with an error of 5% compared to 10% for the DFT algorithm. The most challenging fiber configurations to reconstruct are again very low relative thicknesses and very steep fibers: for d < 0.1 and α < 16◦ , the relative error amounts to hσdi > 1 for both algorithms. For very steep fibers, the error observed in the results of the DFT algorithm increases to up to 80%. The ROFL algorithm achieves a lower error in this case, yet the best case for d = 0.9 still expresses a relative error of hσi ≈ 42%.

#### 3.1.2. Inclination Bias Evaluation

As the orientation vectors were distributed uniformly on a sphere, the frequency of the ground truth inclinations is proportional to the circumference of the cross section of the sphere at the respective height corresponding to a given inclination. In our definition of the inclination angle as α ∈ [− π 2 , π 2 ], the circumference of the cross section of the unit sphere at a respective inclination is proportional to cos(α). In consequence, the inclination frequency p(α) is expected to be proportional to cos(α): p(α) ∝ cos(α).

The inclination histograms obtained from the simulated dataset of uniformly distributed orientation vectors are depicted in **Figure 3**. The histogram of the ground truth inclinations does indeed follow the expected cos(α) proportionality as can be seen from the a · cos(α) curve fit with proportionality factor a in cyan. The proposed ROFL approach achieves a high agreement with the ground truth inclination. On the other hand, the histogram of the inclinations computed by the DFT algorithm reveals a severe lack of in-plane inclinations. Especially orientations with α = 0 ◦ which are the most probable orientations are barely computed at all. Instead of an increasing density toward 0◦ , the histogram displays two symmetric peaks left and right of 0◦ indicated by black arrows which are not present in the ground truth. Another observation lies in the frequency of very steep fibers with respect to the sectioning plane: for |α| > 80◦ , the frequency of DFT inclinations decreases moderately as indicated by blue arrows.

### 3.2. Experimental Data

To demonstrate the working principle of the ROFL algorithm by means of experimental data, one pixel of the stratum sagittale highlighted in **Figure 4A** was chosen. The obtained measurement data from the planar and two tilted measurements after calibration and coregistration are plotted in **Figure 4B**. **Figure 4C** then shows the normalized light intensities and their best fit curves according to the polarimetric model as predicted by ROFL.

The boundaries of the investigated 230 coronal sections are illustrated in the view of the blockface volume shown in **Figure 5A** (Schober et al., 2015; Wiese, 2017). A comparison of the reconstructed histological sections and the corresponding part of the blockface volume is depicted in **Figures 5B,C** utilizing the clipping box view technique implemented in the PLIVIS tool (Schubert et al., 2016, 2018). As the blockface images capture the reflected light during brain sectioning and the transmittance image the transmitted light, they appear inverted to each other. As the time between section mounting and measurement was not constant for all sections, the transmittance differs over the volume as indicated by the white arrows in **Figure 5C**. The reconstruction precision is demonstrated by the smooth reconstruction at the white matter/gray matter transition zones, but also by the fine-grained blood vessel structures highlighted in green in **Figure 5C**.

Volumetric views of 3D-PLI modalities are shown in **Figure 6**. For clarity, the view is the same as in **Figure 5**. The modalities retardation |sin δ| (cf. **Figure 6A**) and relative section thickness d (cf. **Figure 6B**) revealed a strong agreement in most brain regions of the reconstructed volume. Differences, however, were observed in particular for white matter fiber tracts characterized by low retardation values. Those tracts took courses (close to) perpendicular to the sectioning plane according to the fiber orientation map shown in **Figure 6C**. The forceps major (FM) and the stratum sagittale (SS) were worked out as prominent examples for such tracts (cf. **Figure 6**). Ultimately, the obtained fiber orientations (**Figure 6**), in particular in deep white matter structures appeared coherent across neighboring sections. Fiber tracts and pathways were clearly distinguishable from each other and could be traced through the reconstructed volume.

#### 3.2.1. Comparison of ROFL and DFT Algorithms

Major differences between the ROFL algorithm and the DFT algorithm were particularly observed for two cases: low relative section thickness and very steep fibers with respect to the sectioning plane. Therefore, the vector fields of two regions were investigated more closely: the stratum sagittale which is expected to run perpendicular to the sectioning plane and a region at the boundary of white and gray matter. The vector fields underlaid with the retardation maps are shown in **Figure 7**. The plotted two-dimensional, colored lines are representations of the projections of three-dimensional fiber orientation vectors into the respective plane, color-coded by the 3D orientation. The region of interest at the transition zones of white and gray matter shown in **Figure 7A**. Both vector fields are very similar, but the ROFL algorithm results in less inclined orientations than the DFT approach. As the plotted vectors represent the projection of the three-dimensional fiber orientation vectors into the respective plane, longer vectors imply flatter fibers with respect to the respective plane which is the case here. The vector field in **Figure 7B** shows the xz plane perpendicular to the coronal sectioning plane to highlight the robust reconstruction of the stratum sagittale. Minor differences are visible at the boundary of the stratum sagittale but overall both algorithms

FIGURE 4 | Working principle of the ROFL algorithm demonstrated for a single pixel. (A) Transmittance map. The red circle points out the position of the analyzed pixel. (B) Measured light intensities of the planar measurement and tilted measurements to west and east after calibration and registration onto the planar measurement. (C) Normalized light intensities of the planar measurement and tilted measurements to west and east. The dashed lines depict their best fit curves according to the ROFL algorithm. Fit result: <sup>ϕ</sup> <sup>=</sup> <sup>101</sup>◦ , <sup>α</sup> = −60◦ , d = 0.5, R <sup>2</sup> = 0.97.

FIGURE 5 | Reconstructed 3D-PLI volumes. (A) Full blockface volume with black planes representing the boundaries of the analyzed sections shown in (B,C). (B) Partial blockface volume of the histologically analyzed and reconstructed sections. (C) Transmittance volume reconstructed from the histological sections. The green arrows highlight reconstructed blood vessel structures. Note, for the reasons of clarity only vessels with very strong birefringence signals are shown here. The blue arrows indicate white/gray matter transition zones. Note that brightnesses variations pointed out by white arrows occur due to differing times between mounting of the section on the glass slide and the 3D-PLI measurement.

agreed well. As on simulated data the inclination histograms exposed a distinct lacking of in-plane fiber orientations with respect to the sectioning plane for the DFT algorithm, the inclination histograms obtained from experimental data were also investigated. For this purpose one region of interest from the white matter of one section was examined. The fiber orientation

maps (FOMs) resulting from ROFL and DFT are shown in **Figures 8A,B**, respectively. A general glance revealed similar orientation values for both algorithms, which holds true for white and gray matter regions. Due to the considerable amount of noise in the fiber orientations in gray matter, only white matter pixels were considered for the inclination histograms. A detailed analysis of the inclination angle histograms of the white matter regions delineated in **Figure 8C** yielded the following observations (cf. **Figure 8D**): while the frequency of inclination angles drastically decreases from −15◦ to 0◦ and increases again to 20◦ for the DFT algorithm, this lacking of orientations was not observed for the ROFL algorithm. For the regime of high absolute inclination angles, both histograms agree: both have peaks for inclination angles of ≈ ±60◦ (cf. arrows in **Figure 8D**) which result from fibers oriented perpendicular to the sectioning plane in the stratum sagittale (cf. **Figures 8A–C**).

In addition to the potentially biased visual inspection of the vector fields, an unbiased statistical statement about the reconstruction accuracy is given by the mean residual: the mean value of the mean residuals of all pixels is 15% lower for the results of the ROFL algorithm than the results of the DFT algorithm.

#### 3.2.2. Agreement of Model and Data

Since the ROFL algorithm performs a least squares fit in all image pixel, the actual sum of residuals χ 2 (as a result of the optimization process) was accessible. For the DFT algorithm the residual map can also be calculated, yet this requires additional computations which take even longer than the runtime of the algorithm itself. **Figure 9** shows maps of the residuum χ 2 (**Figure 9A**), the relative thickness d (**Figure 9B**) and the fiber orientation map (**Figure 9C**) for a selected region of interest.

As χ 2 is a measure of the difference between the model and the experimental data, it is expected to increase for artifacts such as dust particles which rotate with the polarization filters. Such artifacts were indeed observed in experimental results (cf. **Figure 9A**, where dust particles are indicated by green arrows). For one pixel, a dust particle corrupts one of the 18 obtained values during filter rotation. Still it does not necessarily have a strong effect on the resulting fit values as the corresponding relative thicknesses and fiber orientations do not show signs of artifacts. The χ <sup>2</sup> measure, however, is sensitive to such artifacts.

In addition, we found areas of increased residuals (cf. red arrows in **Figure 9**) we were able to assign to crossing fiber regions. In the fiber orientation map, crossings stood out as the center of an abrupt change of colors (**Figure 9C**). Also, the optimization process resulted in an abrupt change of the relative thickness, even with values of d > 1 for which the relationship between the retardation and inclination is not bijective anymore (**Figure 9B**). Due to the bad agreement between the model and the data, the residuum finally increased strongly (cf. red arrows in **Figure 9**).

### 4. DISCUSSION AND OUTLOOK

We introduced the least-squares algorithm ROFL for the reconstruction of fiber orientation and the extraction of the relative section thickness from measurements with a tiltable specimen stage in 3D-Polarized Light Imaging (3D-PLI). This method requires only one additional assumption besides the polarimetric model: the refractive index of the brain tissue. This represents a substantial improvement as compared to other histological imaging techniques which strongly rely on parameter dependent image processing pipelines to extract three-dimensional fiber orientations. To our knowledge, no other histological imaging technique is currently capable of deriving three-dimensional information from a biophysical model. This compensates for the disadvantage of the difficulties accompanying the 3D re-alignment of serial high-resolution brain section images as required for 3D-PLI anaylsis.

The working principle of the ROFL algorithm was proven for simulated data for which it resulted in a significant improvement of the orientation reconstruction. It was opposed to a previously implemented analytical algorithm based on a discrete Fourier

transformation (DFT). The reconstruction with ROFL became more reliable, in particular for low relative section thicknesses. The noise instability of the DFT algorithm for very flat fibers was already observed by Wiese et al. (2014) and is explainable by the fact that for very flat fibers the gradient <sup>∂</sup> sin <sup>δ</sup> ∂α becomes small (for in-plane fibers with α = 0 ◦ it even becomes zero: ∂ sin δ ∂α <sup>|</sup>α=<sup>0</sup> ◦ = 0). This makes the DFT algorithm prone to noise effects in this case. In fact, the simulations proved that the DFT algorithm is not capable of reconstructing in-plane orientations. On the contrary, the presented fitting algorithm still enabled a reliable reconstruction for in-plane fibers. The random orientations computed by the DFT algorithm for small relative section thicknesses of d < 0.05 were not surprising, as for d < 0.05, the maximal possible retardation value is also |sin δ| = 0.05 and the retardation values of the different tilting positions become almost indistinguishable. In this case, the ROFL algorithm provides a superior reconstruction. In the case of very steep fibers, the ROFL algorithm also outperformed DFT significantly. Still, the reconstruction accuracy decreases strongly for very steep fibers also for the ROFL algorithm. This originates from the signal itself: for α = 90◦ , the amplitude of the sinusoidal signals are almost zero and the measured signals resemble a constant function with random noise (cf. Dohmen et al., 2015 for a theoretical description of this issue). While the reconstruction of out-of-plane fibers remains challenging, orientation vectors with |α| > 80◦ only make up approx. 3% of all possible orientations. Therefore, we can conclude that ROFL accomplishes a very robust orientation reconstruction with an average accuracy of 2◦ for the vast majority of possible fiber orientations for relative section thicknesses of d > 0.2 on synthetic data.

ROFL was applied to 230 consecutive human brain sections subjected to 3D-PLI. The results were 3D reconstructed to demonstrate the robustness and reproducibility of our approach. The obtained modality volumes of retardation, relative section thickness and fiber orientation were characterized by high coherence across the sections. Anatomical structures, such as fiber tracts, white/gray matter borders, vessels could be reconstructed with high precision. The transmittance volume displayed brightness variations, yet the strength of our approach is that the determination of the fiber orientation and relative section thickness is independent of the transmittance due to the normalization. By eyesight, no major differences were observed in the vector fields resulting from the ROFL and DFT approaches. In the stratum sagittale, both approaches resulted in very similar orientations with very minor differences. The computed inclinations were on average α ≈ 65◦ and the computed relative section thicknesses d = 0.6. This combination actually corresponds to the minimum of the orientation reconstruction error of the DFT algorithm obtained from the simulated datasets, which suggests that a strong agreement with the ROFL algorithm can be expected. At the boundary of white and gray matter regions, the ROFL algorithm resulted in less inclined fiber orientations than DFT but also only very minor differences were observable.

The inclination histograms, however, yielded a very important finding: whereas, the frequency of the computed inclination angles decreased significantly for in-plane fibers with α ≈ 0 ◦ for the DFT algorithm, the ROFL algorithm was still capable of reconstructing the whole spectrum of possible inclinations. The same behavior was observed for simulated data, which proves a systematic bias of the DFT algorithm. The ROFL result is also more plausible as for the observed ROI in **Figure 8** the inclination histogram obtained from ROFL agrees far better with the inclination histogram of uniformly distributed orientations. Still, it has to be noted that fiber orientations in a small ROI of the brain cannot expected to be perfectly uniformly distributed due to the convoluted structure of the human brain. The observed inclination differences between ROFL and DFT are hard to observe based on orientation vectors, as the actual inclination differences of up to 10◦ for in-plane orientations can barely be distinguished by eye. As the difference between the measured and predicted light intensities according to the 3D-PLI model also decreases significantly for the parameters obtained from ROFL compared to DFT, ROFL yields more accurate results than DFT also from a purely statistical perspective.

Furthermore, the ROFL algorithm enables a direct way to evaluate the agreement of model and data based on the residuum map which would otherwise have to be computed additionally. This kind of map was exploited to identify measurement artifacts and, even more importantly, crossing fiber regions which were not classifiable by previous 3D-PLI analysis approaches (Axer et al., 2011a; Kleiner et al., 2012; Wiese et al., 2014). Fiber crossings currently pose the greatest complication for 3D-PLI analysis as partial volume effects are still present in brain tissue at voxel sizes of 64 × 64 × 70µm<sup>3</sup> . A voxel containing fibers with different courses (i.e., orientations) will inevitably result in a 3D-PLI measurement composed of superimposed birefringence signals. This might lead to a biased orientation interpretation. To give an example, a voxel comprising crossing out-of-plane and in-plane fibers, the ROFL algorithm will preferably reconstruct the orientation which causes a larger signal, in this case the inplane orientation. To overcome this issue, different strategies are currently being explored: (i) Modeling of crossing fiber structures and subsequent simulations of tilting measurements utilizing the simPLI simulation platform introduced by Dohmen et al. (2015). By this means the behavior of the ROFL algorithm in case of wellknown crossing fiber constellations can be further investigated. (ii) Realizing the idea of oblique imaging for in-plane resolutions at the micrometer scale. Preliminary measurements of a mouse section using a prototypic "tilting" polarizing microscope built up at an optical bench (pixel size: 1.3 × 1.3µm) already revealed the benefit of smaller voxel sizes.

While the obtained parameters show a high coherence across the volume, no quantitative statement about the reliability of the resulting parameters is possible at this point. Hence, future studies need to investigate the uncertainty of the fitted parameters. For Diffusion Tensor Imaging, for example, bootstrapping approaches were used to obtain orientation confidence maps (Jones, 2003; Heim et al., 2004; Whitcher et al., 2008). Especially, the uncertainty of the relative section thickness is of interest, as this parameter is as an indicator for myelin and might even be a reliable measure of the local myelin density. As a matter of fact, the relative section thickness is proportional to the local birefringence 1n which was clearly attributed to the myelin sheath in earlier studies (Goethlin, 1913; Schmidt, 1923; Schmitt and Bear, 1937; de Campos Vidal et al., 1980). Since myelin density information is of outmost interest for the research of degenerative brain diseases characterized by locally altered myelination (e.g., multiple sclerosis), future studies will address the correlation between the relative section thickness and myelin density.

The improved reconstruction comes at the cost of computation time, which increases by a factor of 5,000. Yet, as processing a single brain section takes about 3 min using four compute nodes on JURECA (Jülich Supercomputing Centre, 2016) with the ROFL algorithm in its current implementation, even the computation of whole brains becomes feasible. Still, the computation time could further be reduced by utilizing GPU ressources (Przybylski et al., 2017).

To conclude, the present approach has opened up a new way to determine physical tissue properties from oblique measurements in microscopic 3D-PLI. Cortical, subcortical, and white matter regions could be characterized coherently across brain sections in terms of fiber orientation and birefringence strength. This is a prerequisite for subsequent volume-based connectivity analysis.

#### AUTHOR CONTRIBUTIONS

DS: study design, algorithm implementation, simulation, processing of experimental data, evaluation, writing, discussion.

#### REFERENCES


SM and MS: 3D reconstruction, discussion. NS: 3D visualization. MM: provided the brain sample and anatomical interpretations. KA and TL: writing, discussion. MA: study design, discussion, writing.

#### FUNDING

This project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under Grant Agreement No. 7202070 (HBP SGA1) and Grant Agreement No. 785907 (HBP SGA2).

#### ACKNOWLEDGMENTS

We would like to thank Marcus Cremer and Patrick Nysten for the preparation of the examined brain sections and Felix Matuschke for fruitful discussions about the statistical analysis. We also thank David Gräßel who significantly helped in identifying anatomical structures. The authors gratefully acknowledge the computing time granted by the JARA-HPC Vergabegremium and provided on the JARA-HPC Partition part of the supercomputer JURECA at Forschungszentrum Jülich.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnana. 2018.00075/full#supplementary-material

means of independent component analysis. Neuroimage 49, 1241–1248. doi: 10.1016/j.neuroimage.2009.08.059


Goodman, J. (2000). Statistical Optics. New York, NY: Wiley.


**Conflict of Interest Statement:** 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.

Copyright © 2018 Schmitz, Muenzing, Schober, Schubert, Minnerop, Lippert, Amunts and Axer. 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.

# The Superoanterior Fasciculus (SAF): A Novel White Matter Pathway in the Human Brain?

Szabolcs David<sup>1</sup> \*, Anneriet M. Heemskerk <sup>1</sup> , Francesco Corrivetti <sup>2</sup> , Michel Thiebaut de Schotten<sup>3</sup> , Silvio Sarubbo<sup>4</sup> , Francesco Corsini <sup>4</sup> , Alessandro De Benedictis <sup>5</sup> , Laurent Petit <sup>6</sup> , Max A. Viergever <sup>1</sup> , Derek K. Jones <sup>7</sup> , Emmanuel Mandonnet <sup>2</sup> , Hubertus Axer <sup>8</sup> , John Evans <sup>7</sup> , Tomás Paus ˇ 9,10 and Alexander Leemans <sup>1</sup>

1 Image Sciences Institute, University Medical Center Utrecht and Utrecht University, Utrecht, Netherlands, <sup>2</sup>Department of Neurosurgery, Lariboisière Hospital, APHP, Paris, France, <sup>3</sup>Centre for Neuroimaging Sciences, Institute of Psychiatry, King's College London, London, United Kingdom, <sup>4</sup>Structural and Functional Connectivity Lab Project, Department of Emergency, Division of Neurosurgery, "S. Chiara" Hospital, Azienda Provinciale per i Servizi Sanitari (APSS), Trento, Italy, <sup>5</sup>Department of Neurosciences, Division of Neurosurgery, "Bambino Gesù" Children Hospital, IRCCS, Rome, Italy, <sup>6</sup>Groupe d'Imagerie Neurofonctionnelle (GIN), Institut des Maladies Neurodégératives (IMN)-UMR5293-CNRS, CEA, Université de Bordeaux, Bordeaux, France, <sup>7</sup>Cardiff University Brain Research Imaging Centre, School of Psychology, Cardiff, United Kingdom, <sup>8</sup>Hans Berger Department of Neurology, Jena University Hospital, Friedrich-Schiller University Jena, Jena, Germany, <sup>9</sup>Bloorview Research Institute, Holland Bloorview Kids Rehabilitation Hospital, Toronto, ON, Canada, <sup>10</sup>Departments of Psychology and Psychiatry, University of Toronto, Toronto, ON, Canada

#### Edited by:

Francisco Clasca, Autonomous University of Madrid, Spain

#### Reviewed by:

Miguel Ángel García-Cabezas, Boston University, United States Francesco Sammartino, The Ohio State University, United States

> \*Correspondence: Szabolcs David s.david@umcutrecht.nl

Received: 25 November 2018 Accepted: 07 February 2019 Published: 05 March 2019

#### Citation:

David S, Heemskerk AM, Corrivetti F, Thiebaut de Schotten M, Sarubbo S, Corsini F, De Benedictis A, Petit L, Viergever MA, Jones DK, Mandonnet E, Axer H, Evans J, Paus T and Leemans A (2019) The Superoanterior Fasciculus (SAF): A Novel White Matter Pathway in the Human Brain? Front. Neuroanat. 13:24. doi: 10.3389/fnana.2019.00024 Fiber tractography (FT) using diffusion magnetic resonance imaging (dMRI) is widely used for investigating microstructural properties of white matter (WM) fiber-bundles and for mapping structural connections of the human brain. While studying the architectural configuration of the brain's circuitry with FT is not without controversy, recent progress in acquisition, processing, modeling, analysis, and visualization of dMRI data pushes forward the reliability in reconstructingWM pathways.Despite being aware ofthe well-known pitfalls in analyzing dMRI data and several other limitations of FT discussed in recent literature, we presentthe superoanterior fasciculus (SAF), a novel bilateral fibertract inthe frontal region of the human brain that—to the best of our knowledge—has not been documented. The SAF has a similar shapetothe anterior part ofthe cingulum bundle, but it is locatedmore frontally. To minimize the possibility that these FT findings are based on acquisition or processing artifacts, different dMRI data sets and processing pipelines have been used to describe the SAF. Furthermore,we evaluatedthe configuration oftheSAFwith complementarymethods, such as polarized light imaging (PLI) and human brain dissections. The FT results of the SAF demonstrate a long pathway, consistent across individuals, while the human dissections indicate fiber pathways connectingthe postero-dorsal withthe antero-dorsal cortices ofthe frontal lobe.

Keywords: brain, diffusion MRI, fiber tractography, dissection, polarized light imaging, neuroanatomy, superoanterior fasciculus

### INTRODUCTION

Fiber tractography (FT) based on diffusion magnetic resonance imaging (dMRI; Jeurissen et al., 2017) is widely used for investigating microstructural properties of white matter (WM) fiber-bundles (Alexander et al., 2017), and for mapping structural connections of the human brain (Wakana et al., 2004; Sotiropoulos and Zalesky, 2017). Since the first endeavors of FT

**91**

(Basser et al., 1994, 2000; Conturo et al., 1999; Jones et al., 1999; Mori et al., 1999; Catani et al., 2002), many studies have contributed to the improvement of data quality and the level of detail. Recent advances in MRI hardware and acquisition (Anderson, 2005; Moeller et al., 2010; Setsompop et al., 2013, 2018; Sotiropoulos et al., 2013; Andersson et al., 2016) have improved significantly dMRI data in terms of spatial resolution, angular resolution signal-to-noise ratio, and geometrical distortion reduction. Additionally, notable developments in data processing have reduced errors from the individuals' head motion (Leemans and Jones, 2009), data outliers (Chang et al., 2005; Pannek et al., 2012; Collier et al., 2015; Tax et al., 2015), eddy currents (Andersson and Sotiropoulos, 2015, 2016; Andersson et al., 2016), EPI distortions (Andersson et al., 2003), Gibbs-ringing (Perrone et al., 2015; Kellner et al., 2016; Veraart et al., 2016), and other data artifacts, like drift of the diffusion signal (Vos et al., 2017). Besides the traditional diffusion tensor model, more sophisticated methods have been developed to resolve crossing fibers (Lin et al., 2003; Tuch, 2004; Tournier et al., 2004; Wu and Alexander, 2007; Dell'Acqua et al., 2013; Tax et al., 2014; Jensen et al., 2016). Together, all these improvements allow for more reliable FT results, thereby resolving complex fiber architecture, smaller branches of fiber bundles or minor fiber-pathways (Tournier et al., 2008). Although dMRI based tractography is a promising technique, there are many well-known pitfalls and limitations in acquiring and analyzing dMRI data (Jones and Cercignani, 2010; Jones et al., 2013b; O'Donnell and Pasternak, 2015). Thus, dMRI is an indirect approach for measuring the underlying WM pathway properties, as it is an approximation of the average diffusion within a given voxel. Modeling the diffusion results in a proxy of main directions, and thus the exact architecture of the pathways cannot be determined unambiguously. The modeling errors can cause tracts to stop prematurely or jump from one WM structure to another, resulting in false negative or false positive connections. The strength of connectivity, when assessed via probabilistic FT, is unreliable due to its sensitivity to data quality (Mesri et al., 2016).

Despite recent efforts to increase the accuracy of FT by either pruning the tractograms using anatomical information (Smith et al., 2012; Roine et al., 2015) or by including microstructural information to disambiguate between pathways (Smith et al., 2013), the International Society for Magnetic Resonance in Medicine (ISMRM) Tractography Challenge in 2015 demonstrated that some data-processing pipelines could result in large errors as the reconstructed tracts produced an average ratio of 1:4 in false positivefalse negative connections and 45% in bundle overlap when compared with pre-defined, ground-truth streamlines (Maier-Hein et al., 2017).

In another study, the sensitivity (true connections) and specificity (avoidance of false connections) of dMRI tractography was investigated by combining tracer studies and high quality dMRI data (Thomas et al., 2014). This research showed a sensitivity-specificity trade-off, i.e., the number of true connections increased simultaneously with the number of false connections. Additionally, the anatomical accuracy of the tracts depended on the studied pathway, the acquisition and the processing pipeline. This stresses the fact that comparing results across studies can suffer from the differences in acquisition settings and processing methods. Therefore, reconstructing fiber pathways from dMRI FT should be performed with great care and additional support from other methods is welcome.

Nonetheless, being aware of the recognized issues, we present the description of a novel fiber tract, the superoanterior fasciculus (SAF), in the frontal lobe which—to the best of our knowledge—has not been documented before with dMRI based tractography with this level of detail. The tract is slightly curved and follows the same arc profile of the cingulum bundle, but is located more frontally and above the frontal part of the cingulum; it can be found in both hemispheres.

Traditionally, the description of new WM fiber bundles went through several stages. For example, WM histology and bundle degeneration studies (Dejerine and Dejerine-Klumpke, 1895) and dissection (Ludwig and Klingler, 1956; Heimer, 1983) studies supported the presence of a structure in post mortem brains, which was confirmed later via dMRI-based FT in the living human brain. Conturo et al. (1999) laid the foundations for selection and analysis of tracts based on ''region of interest (ROI)'', which opened up a new era of tract-based investigations from neuroscience (Lebel et al., 2008; Thiebaut de Schotten et al., 2011a, 2012) to clinical applications (Thiebaut de Schotten et al., 2005; Bartolomeo et al., 2007; Deprez et al., 2012). The access to large samples made the construction of in vivo WM tract atlases and guidelines feasible (Mori et al., 2005; Wakana et al., 2007; Hua et al., 2008). In this work, we reconstruct the proposed structure using specific ROIs along with modeling of the diffusion MR signal based on constrained spherical deconvolution (CSD; Tournier et al., 2004; Tournier et al., 2007).

To minimize the possibility that these findings are based on acquisition or processing artifacts, we used datasets from different projects: young Adult human connectome project (HCP; Glasser et al., 2013; Sotiropoulos et al., 2013; Van Essen et al., 2013), Avon Longitudinal Study of Parents and Children (ALSPAC; Golding, 2004), Multiple Acquisitions for Standardization of Structural Imaging Validation and Evaluation (MASSIVE; Froeling et al., 2017), an in-house dataset (Jeurissen et al., 2011), and robust processing pipelines to show these trajectories. While the main focus of this research was based on dMRI, we complemented our findings with other non-MRI techniques, such as polarized light imaging (PLI) microscopy (Axer et al., 2011) and human brain dissections. Preliminary results of this work on the SAF have been presented at the 2015 ISMRM meeting in Toronto, ON, Canada (Heemskerk et al., 2015).

#### MATERIALS AND METHODS

In the methods section, we describe first the datasets that were used for the investigation. We then explain the dMRI processing, FT and analysis steps. Finally, we describe the PLI and dissection methods.

### Diffusion MRI

#### Datasets

Four dMRI datasets acquired on different platforms were used to study the tract. The datasets represent a large diversity in data quality and number of samples, which helps to generalize the conclusions.

#### **Dataset 1**

Minimally processed dMRI data from the HCP S500 release were used. Briefly, the data were acquired with a 1.25 mm isotropic voxel size; three shells with b = 1,000, 2,000 and 3,000 s/mm<sup>2</sup> in 90 DW volumes and six non-weighted images per shell. Every participant for whom all the 90 b = 1,000 s/mm<sup>2</sup> and 90 b = 3,000 s/mm<sup>2</sup> images were available along with all 18 b = 0 s/mm<sup>2</sup> volumes, and were not listed among participants with known anatomical anomalies or data quality issues were included in the analysis. The IDs of the excluded participants are listed on the HCP wiki page (HCP wiki, 2018). The selection resulted in 409 healthy participants (243 females and 166 males), between the age of 22 and 36 years. Low b-value images were analyzed with DTI and CSD models separately using 90 b = 1,000 s/mm<sup>2</sup> and nine unweighted volumes, while the higher b-shell was processed with CSD only using 90 b = 3,000 s/mm<sup>2</sup> and 18 unweighted volumes. For a subset of 10 participants we also analyzed the b = 1,000 s/mm<sup>2</sup> , 2,000 s/mm<sup>2</sup> and 3,000 s/mm<sup>2</sup> shells with both CSD and DTI models.

#### **Dataset 2**

Ten dMRI datasets were used from the ALSPAC cohort, each consisting of 30 diffusion gradients with b = 1,200 s/mm<sup>2</sup> , and 2.4 mm isotropic voxel size as described by Björnholm et al. (2017).

#### **Dataset 3**

The dataset is described in the work of Jeurissen et al. (2011). In summary, one participant's dMRI consisted of 60 diffusion directions, 2.4 mm isotropic voxel size (resampled to 1 mm isotropic) and b = 3,000 s/mm<sup>2</sup> .

#### **Dataset 4**

The dataset from the MASSIVE (Froeling et al., 2017) acquisition was used. In short, we used a subset of 22 b = 0 s/mm<sup>2</sup> volumes and 250 b = 3,000 s/mm<sup>2</sup> volumes and an isotropic resolution of 2.5 mm<sup>3</sup> of one participants.

#### Image Processing and Signal Modeling

All datasets were preprocessed with ExploreDTI v.4.8.6 (Leemans and Jones, 2009; Leemans et al., 2009), except for the HCP dataset, which was already preprocessed (see Dataset 1 below). Further modeling (DTI and CSD) and defining ROI configurations were also performed with ExploreDTI.

#### **Dataset 1**

Briefly, FSL (Jenkinson et al., 2012) tools topup (Andersson et al., 2003) and eddy (Andersson and Sotiropoulos, 2016) were used to correct for head motion and geometrical distortions arise from eddy currents and susceptibility artifacts. Finally, the DWIs were aligned to the structural T1 image. For the full description, see the work of Sotiropoulos et al. (2013). DTI estimation of the low b-value shell for the full population was performed using REKINDLE (Tax et al., 2015). The fiber orientation distribution (FOD) in each voxel was estimated using CSD (Tournier et al., 2004, 2007) with the recursive calibration method (peak ratio threshold = 0.01; Tax et al., 2014) with maximum harmonic degree of Lmax = 8, in the high and low b-value shells separately. For a subset of 10 subjects, we performed additional analyses to test the effects of b-value and choice of diffusion modeling. For this, we used the following settings: (a) CSD with b = 2,000 s/mm<sup>2</sup> ; (b) CSD with b = 1,000 s/mm<sup>2</sup> ; (c) CSD with Lmax = 6; (d) CSD with calibration of the response function based on FA > 0.8; and (e) DTI estimation using REKINDLE.

#### **Dataset 2**

Motion-distortion correction and Gaussian anisotropic smoothing (Van Hecke et al., 2010) were performed with ExploreDTI. The FOD in each voxel was estimated using CSD with recursive calibration (peak ratio threshold = 0.01) and Lmax = 8.

#### **Dataset 3**

Motion-distortion correction was performed with ExploreDTI. The dataset was resampled to 1 mm isotropic voxel size to increase the level of detail (Dyrby et al., 2014). FODs were estimated using CSD with recursive calibration (peak ratio threshold = 0.01), Lmax = 8.

#### **Dataset 4**

Dataset 4 was motion, distortion corrected and resampled to 1 mm isotropic voxel size. Similar to the previous datasets, the FOD was estimated using CSD with recursive calibration (peak ratio threshold = 0.01), Lmax = 8.

#### Fiber Tractography

For all four datasets, the deterministic FT framework of Jeurissen et al. (2011) was used for the CSD approach, with the following parameter settings: FOD threshold = 0.1, angle deviation = 45◦ , step size = 1 mm and minimal tract length = 20 mm. Wholebrain FT was performed with uniform distribution of seed points defined at a 2 mm resolution. For DTI based FT, the settings are: step size = 1 mm, minimal tract length = 20 mm, angle deviation = 45◦ , FA seed and tracking threshold of 0.2 and uniform seed point resolution of 2 mm.

#### ROI Configuration for Reconstructing the SAF

Automated, large-scale tract selection on the HCP dataset was obtained by atlas-based tractography segmentation (Lebel et al., 2008). On the template dataset, we defined 4 Boolean ''AND'' and ''NOT'' ROIs: 2 axial AND, 1 sagittal NOT and 1 coronal NOT ROIs (**Figure 1**). The 2 axial AND ROIs were placed as follows: the first one was located at the height of approximately half the genu of the corpus callosum (CC) and the second ROI was placed 10 slices (equals to 12.5 mm in the HCP dataset) superior. Both AND ROIs included the medial frontal area and excluded the cingulum. A NOT ROI was placed midsagittal to exclude fiber pathways from the CC. Additionally, one coronal NOT ROI was placed posterior of the frontal lobe and below the CC to exclude fibers from the inferior fronto-occipital fasciculus (iFOF).

FIGURE 1 | Positions of the two AND (in green) and two NOT (in red) regions of interest (ROIs) are shown with a direction encoded color (DEC)-FA map. This ROI configuration was used for the tract selection of the superoanterior fasciculus (SAF) in both hemispheres.

#### Consistency Evaluation of the SAF **Within Cohort**

For each of the 409 participants from the HCP cohort, the visitation mask of the reconstructed SAF was normalized to the Montreal Neurological Institute (MNI) template (Fonov et al., 2011), using the transform functions provided by the HCP team. The resulting map was thresholded at 1 and binarized. The group map displays the fraction of participants for whom the tract mask of the SAF is present in each voxel of the MNI space.

#### **Across Acquisition Protocols**

The tractography results of participants from the four different datasets were compared to investigate the consistency of the SAF configuration across different acquisition methods.

#### **Across Processing Settings**

We evaluated the effects of different processing settings on the tractography results for 10 subjects of the HCP dataset.

#### Polarized Light Microscopy

PLI is a microscopy method that uses the birefringent properties of myelin sheath to quantify the main fiber orientation in histological sections (Larsen et al., 2007). In this work, we used data previously involved in the investigation of the anterior cingulum (Axer et al., 2011). One of the studied brains was sectioned as to expose the area of interest but was not investigated previously. For details on acquisition and processing see the corresponding study (Axer et al., 2011). Briefly, a 4% aqueous formalin-fixed human brain was macroscopically dissected and a 1.5 cm thick slab of the medio-frontal brain including the anterior cingulum bundle was cut in four blocks. Each of these four blocks were serially sliced and analyzed. Sections were obtained with a cryostat microtome (CM3050 S, Leica Microsystems, Bensheim, Germany) at a thickness of 100 µm. Aqueous mounting medium AquatexTM (Merck, Darmstadt, Germany) was used for mounting. The histological sections were placed between two coupled crossed polars which can be rotated. Birefringence in the tissue is able to twist some of the light so that it can pass through the second polarizer and be imaged. The orientation of the nerve fibers influences the transmission of plane-polarized light at different azimuths. A CCD camera (Axiocam HR, Carl Zeiss, Göttingen, Germany, basic resolution of 1,388 × 1,040 pixel) was used to capture the light. The brain slices were digitized under azimuths from 0 to 80◦ using two polars only. These sequences were used to estimate the inclination of fibers (in z-direction). The same slices were digitized under azimuths from 0 to 160◦ in steps of 20◦ using a quarter wave plate additionally. These sequences were used to estimate the direction of the fibers in xy-direction.

### Dissection (Paris)

Post-mortem brain dissection was performed in the laboratory of neuroanatomy of Lariboisière Hospital in Paris, France by two neurosurgeons (EM, FCorr).

Five human cerebral hemispheres, obtained from fresh autopsy, were fixed in 10% formalin solution for at least 3 weeks. These specimens were frozen at −18◦C for 2 weeks and then unfrozen at room temperature. This process (Klingler technique's; Klingler, 1935) induces water crystallization in the brain tissue that, by spreading along the WM fibers, facilitates visualization and dissection of the subcortical WM fibers.

The dissection was performed with light microscopy (OPMI pico, Carl Zeiss, Inc., Oberkochen, Germany) and video recorded with a high-resolution digital camera (Karl Storz GmbH, Tuttlingen, Germany).

Initial observation of the configurational shape of sulci and gyri of the mesial cortical surface of frontal lobe was always performed before dissection. The dissection always started from the cingulate fibers extending to the entire mesial and lateral surface of the frontal lobe.

#### Dissection (Trento)

Dissections below were performed by three neurosurgeons (SS, FCors, and AB) in the context of the Structural and Functional Connectivity Lab (SFC-Lab) Project, which was approved by the Ethical Committee of the Azienda Provinciale per i Servizi Sanitari (APSS) of Trento, Italy.

Three right hemispheres were prepared for Klingler's dissection (Klingler, 1935) according to the protocol previously reported (Sarubbo et al., 2015; De Benedictis et al., 2014, 2016). It starts with an injection of formalin 10% in carotids and vertebral arteries, then an immersion in formalin 10% for 40 days, and finally the progressive freezing at −80◦ and de-freezing procedure. After the first de-frost process and removal of arachnoids and vessels, we started the micro-dissection under microscope with 4× magnification with wooden spatulas, approaching at the medial and ventral WM of the frontal lobe from the latero-dorsal cortical surface and leaving intact the gray matter at the tip of the gyri exposed to highlight the territories of terminations.

### RESULTS

### Tract Description

The architectural configuration of the SAF is shown in **Figure 2** where dataset 3 was used with CSD based FT. The complex fiber architecture in the frontal area can be appreciated from the FODs overlaid on the sagittal view in **Figures 2A–C** with the dotted yellow line identifying the interface between regions with locally different dominant fiber populations [''blue'' vs. ''green'' on the principal direction encoded color (DEC) map]. A medial and frontal view of the bilateral tract configuration is given in **Figures 2D,E**, respectively. Note that the SAF resembles the trajectory of the cingulum bundle (in red) but it is located more anteriorly, superiorly, and laterally.

The SAF group composite map of the HCP participants (dataset 1) using CSD modeling with b = 3,000 s/mm<sup>2</sup> is shown in **Figure 3**. Again, the bilateral tracts follow a similar trajectory as the cingulum, but more frontally [i.e., in front of the cingulate sulcus or within the superior frontal gyrus (SFG)] and slightly more laterally. The SAF pathways appear to spread from the rostrum of the CC to the ascending ramus of the cingulate sulcus and are medial

FIGURE 2 | Location of the SAF: the axial (A) and sagittal (B) view of part of the region that intersects the SAF (blue region indicated by the arrows). The complex fiber architecture in this area can be appreciated from the fiber orientation distributions (FODs) overlaid on the sagittal view in (C) with the dotted yellow line identifying the interface between regions with different dominant fiber populations ("blue" vs. "green" on the DEC-FA map). Sagittal (D) and coronal (E) views of the SAF configuration. For displaying purposes, only one of the bilateral SAFs is shown. The cingulum bundle is shown in red to provide anatomical reference. Grayscale background is the T1-weighted magnetic resonance imaging (MRI) of the subject.

fraction of subjects for which the tract is present in a given voxel.

to the corona radiata. **Supplementary Figures S1, S2** show the population level map using the same ROIs, but with the 1,000 s/mm<sup>2</sup> b-value dataset and modeled with DTI and CSD, respectively.

#### Tract Consistency

#### Across Participants

The percentage overlap map in **Figure 3** and **Supplementary Figure S2** show that the main portion of the SAF is present for 90% of the participants, which is in the same range as other tracts (Thiebaut de Schotten et al., 2011b). The extent of the reconstructed SAF varied substantially across participants (see **Figure 4**). For some, a very extended and broad structure was found (**Figure 4**, bottom row), while in others we found only a narrow portion of the SAF (**Figure 4**, top row). The length of the streamlines also varied across participants, but also within the same individual. It is often observed that streamlines stop abruptly due to the absence of a corresponding FOD peak larger than the predefined threshold.

#### Across Data Acquisitions

**Figure 5** reveals that the SAF was present in all samples that we analyzed. This observation excludes the possibility that the SAF is merely a scanner-specific artifact. While there are differences in acquisition settings between the different acquisition sites, there is no impact on the extent of the SAF. An example is shown in **Figure 6**, where the effect of using different b-values for the data acquisition on the SAF reconstruction is depicted. Low b-values result in less complex FODs leading to shorter and fewer streamlines, but CSD based modeling can still reveal the SAF configuration. CSD based tract coverage from low b-value data shows large overlap with the high b-value based CSD tracts, revealing that proper modeling of crossing fibers is imperative for reconstruction of the SAF. This can be appreciated as well in the population map in **Supplementary Figure S2**, in which DTI based FT cannot deal with complex fiber pathway configurations.

#### Across Processing Strategies

Modeling and processing steps also affected the extent of the reconstructed SAF (**Figure 7**). DTI analysis (**Figure 7A**) resulted in either no pathways at all or in some spurious ones. CSD-based modeling showed either shorter SAF pathways (Lmax = 6, **Figure 7C**) or minor variations (FA calibration, **Figure 7D**). According to the population-level maps in **Supplementary Figure S1**, DTI-based reconstructions of the SAF cover only the most frontal parts of the tract, which is also depicted in **Figure 2** with the blue color. In this location of the SAF, the locally dominant direction of diffusion is tangential with the SAF, irrespective of any other crossing pathways. Clearly the biggest impact on showing the SAF configuration was the modeling of crossing fibers, in our example using CSD over DTI.

#### PLI

**Figure 8** shows the ensemble of the four blocks color-coded by the main in-plane orientation per voxel. In the ROI of where the SAF is expected, a variety of colors and therefore locally dominant orientations were present. Nonetheless, a similar shape (indicated by the white arrows) as with tractography can be appreciated.

#### Dissection (Paris)

We were able to isolate, in four of the five specimens, a series of fibers spreading over from the cingulate fibers that seemed to have the same anteroposterior orientation as the

FIGURE 4 | Consistency within the HCP cohort. The right SAF and cingulum is depicted in sagittal view for nine subjects from the HCP dataset highlighting the large variability in extent of the SAF. The cingulum is shown in red to provide anatomical reference. Tracts are plotted on top of the DEC-FA map.

FT results of the SAF. These fibers were located above the cingulum in correspondence with the frontal gyrus (**Figure 9**). The careful dissection of the fiber complex revealed several different fibers: U-shaped fibers from cingulate, U-shaped fibers belonging to SFG and more lateral fibers from corona radiata (**Figure 9**). All of them had a vertical orientation; therefore, no correspondence to the anteroposterior orientation of the putative SAF.

### Dissection (Trento)

After exploration of the sulco-gyral anatomy (**Figure 10A**) we proceeded with a gentle peel out of the gray matter of the depth and lateral surface of all the gyri of the dorsal portion of the frontal and parietal lobes (**Figure 10B**). Then we removed the U-fibers connecting the most anterior thirds of the middle and superior frontal gyri, at the level of the superior frontal sulcus (SFS; **Figure 10C**). After a progressive and cautious removal of the U-fibers connecting the lateral and medial cortices facing at the border of the SFS, we exposed a thin layer of fibers with an antero-posterior course (in all the three hemispheres dissected; **Figure 10D**). We followed back with a gentle microdissection this layer of fibers exposing the posterior territories of terminations within the dorsal third of the pre-central gyrus. Anteriorly we exposed this thin bundle with an arching course, that follows the physiologic curvature of the frontal anterior cortices, and we exposed the territories of terminations within: the frontal pole, the fronto-orbital cortex and the fronto-orbitolateral cortex.

FIGURE 5 | Consistency across acquisition protocols. For all four datasets an example of the SAF configuration is shown in sagittal view: right SAF and right cingulum for (A) dataset 1 and (B) dataset 3. Left SAF and left cingulum is shown for (C) dataset 2 and (D) dataset 4. The cingulum is shown in red to provide anatomical reference. Tracts are plotted on top of the DEC-FA map.

It was not possible to follow these fibers further back (i.e., behind the central sulcus) in any of the specimen because: (1) we did not find clear signs of continuity; and (2) the dense and intermingled crossing area with the vertical fibers ascending and descending to the pre- and post-central gyri.

This thin layer of fibers is located in the most medial and dorsal portion of frontal WM and connects the postero-dorsal with the antero-dorsal cortices of the frontal lobe (**Figure 11A**). We performed also the microdissection of the most ventral components of the superior longitudinal fascicle (namely, SLF II and III; **Figures 11B–D**), accordingly to the recent description by Fernández-Miranda et al. (2015); these SLF components occupy the mediodorsal and ventral WM of the frontal lobe and are clearly

FIGURE 7 | Consistency of the SAF configuration across different modeling and processing strategies for the same subject of the HCP cohort. In (A), the DTI model was used, resulting in only few and likely spurious tract pathways, whereas the CSD models in (B–D) produce similar and plausible SAF trajectories, but with differences due to processing settings: (B) recursive calibration and Lmax = 8; (C) recursive calibration and Lmax = 6; and (D) FA calibration and Lmax = 8. Tracts are plotted on top of the DEC-FA map and shown sagittal.

distinguishable from the SAF fibers by location, direction and course.

### DISCUSSION

By taking advantage of advanced dMRI methodology, we have identified a consistent bilateral pathway, referred to as the SAF, via FT in the frontal lobe, which has not been described before (Catani et al., 2012; Makris et al., 2005; Thiebaut de Schotten et al., 2011b). Reproducibility across multiple participants, different data samples and acquisition settings boosted our confidence that this finding is not based on imaging artifacts.

FT shows that the SAF can be consistently reconstructed with FT. However, the presence of a long association fiber is not fully supported by dissection, which necessitates hypothesizing about why this pathway is still found with FT. An explanation could be that a series of consecutive U-shaped fibers that connect adjacent gyri emerge to form a long pathway (see **Figures 8**, **12**). This possibility is supported by the work of Maldonado et al. (2012), where they propose that the dorsal component of SLF is primarily composed of U-shaped fibers. However, dissection also revealed a distinct WM structure with the same orientation as the proposed fiber bundle after the explicit removal of U-fibers. One may argue that the fibers are part of the dorsal component of the SLF system (i.e., SLF I), but considering the common anatomical definitions and the territories of termination of the SLF, it is unlikely the case. **Figures 11C,D** highlights the different and well-known components of the SLF. Important to note that this regions is challenging for FT and dissection as well due to the high number of crossings with the descending and ascending (i.e., vertical) pathways.

Independent investigations are still necessary to decide whether the suggested structure is indeed a new one or a component of an already defined tract. A possible description based on new nomenclature that it is part of the large dorsal system of horizontal fibers connecting frontal and parietal cortices as a branch of the mesial longitudinal system (MesLS; Mandonnet et al., 2018), which can be further divided into

FIGURE 8 | Polarized light imaging (PLI) results. (A) Four separate blocks of the same brain sliced in sagittal plane, where the fiber orientation maps are depicted correspond to the (B) color scheme circle. White arrows are indicating the location of the proposed pathway. Part (C) shows an example of formalin fixed human cadaver brain.

an inner branch, which is the cingulum per se and an outer branch, the SAF. Note, such an outer branch is close to what Wang et al. (2016) named as supracingulate or paracingulate pathways. **Supplementary Figure S3** shows the population of the SAF map CSD based modeling and high b-value data along with the cingulum from the JHU-ICBM DTI based WM atlas (Mori et al., 2005; Wakana et al., 2007; Hua et al., 2008) demonstrating that the SAF is not part of the cingulum.

While the agreement between dissection and FT is imperfect, it is important to describe the structure and recognize its occurrence. First, as this is a consistent finding, other researchers will likely discover and investigate the presence of this FT structure. Second, investigating the structural properties of the sub-compartments can improve the understanding of the U-fibers and rule out false structures, which were accidentally reconstructed from U-fibers. Third, this adds to the discussion of the value of FT, thereby raising additional awareness of the pitfalls of FT (Jones and Cercignani, 2010; Parker et al., 2013; O'Donnell and Pasternak, 2015), and the concerns in the field of connectomics (Hagmann et al., 2010; Fornito et al., 2013; Zalesky et al., 2016).

In the following sections, we will discuss that our FT results are plausible from the point of dMRI-based tractography, place our results in the context of other FT and dissection studies, provide suggestions for future investigations, and stimulate the debate about the validity of dMRI-based FT.

### Diffusion MRI

Recent developments in acquisition and data processing boosted the reliability of dMRI and increased the inherently low accuracy of mapping WM pathways with FT. The investigations in solving crossing fibers were imperative, as the reported structure would not or hardly be detected since the main diffusion direction here is along the forceps (see **Figure 2**). However, one of the pitfalls of dMRI is that it is an indirect measure of the underlying WM pathways and reflects the net displacement of water molecules along all structures within a large voxel. Different structural architectures can lead to the same diffusion profile (Jbabdi and Johansen-Berg, 2011), making it hard to reconstruct unambiguously the intrinsic fiber arrangement.

#### Impact of Modeling

The choice of diffusion signal modeling was critical as the SAF cannot be revealed with DTI. The tract is formed mostly from locally non-dominant fiber populations. Therefore, it can be mapped only with techniques, which can resolve crossing fibers. The single location where the tract is locally dominant is the most anterior part with up-down orientation, as depicted in blue in **Figures 1A–C**, and is also the only location revealed by DTI based population level tract mapping shown in **Supplementary Figure S1**. Following this line of reasoning, yet a system of non-dominant, interwoven fiber bundles are to be found with the help of high data quality and advanced modeling.

Another example of mapping a challenging fiber population without any well-known underlying architecture is showcased in the MyConnectome project (Poldrack et al., 2015), featuring longitudinal, high quality MRI data of a single adult male. In a small region of the CC, an anterior–posterior oriented set of fibers are dominating locally, yet the left-right oriented CC fibers are still present and continuous. This rare finding is not expected in the CC, but is remarkably strong in that

FIGURE 9 | Dissection results in two different hemispheres and color overlays representing the proposed structure. Top row: (A) the original image and (B) possible trajectory of the SAF overlaid in blue, the cingulum is indicated in red as the anatomical reference. Bottom row: (C) original image and (D) several U-shaped fibers highlighted in red, which may form part of the SAF via FT.

the front-back oriented diffusion was larger than the left-right orientated diffusion and, hence, changed the color in the DEC map from red to green. Similar to the SAF, it is difficult to investigate parts of pathways where the underlying FOD profiles exhibit non-dominant peaks tangential to the pathway of interest.

#### Tract Consistency

CSD tractography results obtained with low and high b-value data show a high similarity of 90% for the core of the SAF, which is outstanding given the fact that WM fiber bundles vary in their size and position, especially for this is a small and inferior bundle. Compared with other larger bundles, the similarity is in line with previously reported overlap values of >90% (cingulum), >90% (core CC), >75% (CC), and >75% (iFOF; Thiebaut de Schotten et al., 2011b).

The presence of the SAF across different data samples and acquisition settings is indicative that the results are not merely based on acquisition artifacts. As FODs are sharper, more and longer streamlines are found at higher b-values as this will aid the resolving power of the different fiber populations within a voxel. Although we can show the SAF in case of all the participants, finding the complete trajectory is not trivial. The SAF can be incomplete, for example, when the FOD only contains the main peak of the crossing fibers and does not contain the minor second peak. While we were able to determine the main part of the SAF configuration, the dorsal terminations were hard to obtain, which was also an issue during dissection. Future studies, especially at higher spatial and angular resolution, may alleviate this question.

#### Validation

A conceptual limitation of FT (Jeurissen et al., 2017) is the inability to resolve the functionality and directionality of anatomical links, where the existence of a physical connection is necessary to complete the cortical functions. Such functional confirmation can arise from the use cortico-cortical (Matsumoto et al., 2004, 2006, 2012) or axono-cortical evoked potentials (Yamao et al., 2014; Mandonnet et al., 2016), generally obtained with invasive techniques involving implants or electrodes during surgery.

The PLI results show an overall pathway that can resemble the presented structure; however, it is not clearly one consistent pathway. A drawback of PLI, however, is that it cannot distinguish between differently oriented crossing fibers, making it hard to be unambiguous about the exact pathway.

On the one hand, the brain dissections in this work could not clearly verify the existence of the SAF. On the other hand, the dissection results could not rule out the existence of the SAF given the limitations of this approach, especially in regions where multiple fiber systems exhibit crossing configurations. The brain data used for the FT results in this work are obtained from

FIGURE 10 | (A) The main sulci and gyri of the lateral frontal region of a right hemisphere are indicated, after removal of vessels, pia mater, and arachnoid membranes. (B) Decortication of the SFG and MFG allows identifying the U fibers under the SFS (blue pin, pink circle). (C) Exposition of U-fibers between di SFG and MFG (pink arrows). (D) Exposition of a thin layer of fibers with an antero-posterior course under the SFS (pink arrow). CS, central sulcus; IFG, inferior frontal gyrus; IFS, inferior frontal sulcus; MFG, middle frontal gyrus; Op, opercular part of the IFG; Orb, orbital part of the IFG; PreCG, precentral gyrus; PreCS, precentral sulcus; PostCG, postcentral gyrus; S, superior; SFS, superior frontal sulcus; SFG, superior frontal gyrus; Tr, triangular part of the IFG.

healthy volunteers. Hence, ex vivo investigations (dissection and PLI) on the same brains were not feasible. As such, a direct comparison was not possible.

### Complexity of Fiber Bundles

Interaction between FT and brain dissection methods drives both the validation of tracts, but also the discovery of true fiber bundles, thereby increasing our understanding of the design of the brain's architecture (Yeatman et al., 2014; Meola et al., 2015; De Benedictis et al., 2016; Wu et al., 2016b; Hau et al., 2017). Over the years, more and more (parts of) WM bundles were revealed using FT, which were previously undetectable with the formerly existing techniques (Jbabdi and Johansen-Berg, 2011) such as the lateral projections of the CC, IFOF or the Aslant fiber bundle (Thiebaut de Schotten et al., 2012). However, the validity of some of these tracts is still debated. Yeatman et al. (2014) used tractography to rediscover the vertical occipital fasciculus (VOF); a bundle that caused controversies among neuroanatomists in the 19th century. These tractography findings were further confirmed via dissection by Wu et al. (2016c). In this context, the identification of the SAF with FT raises similar concerns and fuels the ongoing debate about validation with other approaches.

Other researchers have also proposed new or redefined WM structures. Track density imaging (TDI) on 7T dMRI (Calamante et al., 2010, 2011, 2012) recently demonstrated finer details of thalamocortical connections (Choi et al., 2018) and revealed the fiber system of the septum pellucidum area (Cho et al., 2015). In both studies the super-resolution ability of the voxelwise fiber count was used to generate images with high anatomical contrast and therefore exposed the aforementioned bundles. However, due to the noise sensitivity of the technique, the validity of structures revealed by TDI remains an open question (Dhollander et al., 2012, 2014).

In addition to new bundles, several studies have shown subcompartments of known WM pathways, which are validated by histology. These multi-component bundles such as the uncinate fasciculus (Hau et al., 2017), iFOF (Sarubbo et al., 2013; Caverzasi et al., 2014; Hau et al., 2016; Wu et al., 2016a) and SLF (Makris et al., 2005; Kamali et al., 2014) consist of a complexity of pathways that together form a united bundle. De Benedictis et al.

FIGURE 11 | Dissection of a right hemisphere according to Klingler's technique. (A) Exposition of the most ventral components of the anterior indirect component of the SLF (SLF II and SLF III; green area). (B) Thin layer of fibers connecting the postero-dorsal with the antero-dorsal cortices of the frontal lobe (white area, red arrow). (C) Dissection of the main components the SLF: indirect anterior (green area), indirect posterior (blue area); direct (yellow area). (D) U-fibers (blue pins), dorsal frontal fibers (white pins), frontal terminations of the AF, respectively in the ventral premotor cortex (green pins) and the pars opercularis of the IFG (red pins). AF, arcuate fasciculus; SLF, superior longitudinal fascicle.

(2016) used a microdissection approach to reveal and validate the presence of both homotopic as well as heterotopic fibers in the anterior half of the CC.

The cingulum, which has a similar shape and is in close approximation of the SAF, is also a multi-component bundle. The cingulum bundle consists of many short fibers as well as longer fibers that together have many different connections (Bajada et al., 2017). The cingulum is usually depicted as a continuous structure using FT although thorough research is showing a division in at least three subparts, each with their own distinct diffusion metrics (Jones et al., 2013a; Wu et al., 2016b). Given the similarity with the cingulum and the knowledge that many bundles have complex and multicomponent fiber organization, it seems plausible that the SAF also has a multicomponent organization.

#### Future Directions

This study describes a new pathway identified with dMRI based FT, called the SAF, which may have several implications in in vivo neuroimaging studies about structural brain connectivity. If the SAF is considered to be a genuine new brain pathway, further research is needed to better understand the origin and the function of this structure. Investigation of this area at higher spatial resolution could be one avenue as recent work has shown that this coincides with a better understanding of complex crossing fiber structures (Schilling et al., 2017). This approach would then ideally be combined with dissection. We can expect more, previously not shown and locally non-dominant structures to be discovered as Jeurissen et al. (2013) showed that most of the WM in the human brain contains multiple populations. Other imaging techniques, such as optical coherence tomography (Huang et al., 1991) or optogenetics (Deisseroth, 2010) may provide complementary information that could verify the existence of the SAF.

We observed a large variability in the cross-sectional area of the SAF and future analysis of the structure should also entail examining common properties related to microstructure, shape, demographics, and even changes in pathological conditions. Furthermore, it is well-known that the frontal lobe is generally a challenging region to investigate by most MRI methods, because of the susceptibility induced distortions, which are further emphasized by EPI sequences caused by the air-tissue interfaces of the sinuses.

In the past, we have seen examples in which part of anatomical knowledge is abandoned from the mainstream, e.g., in case of the VOF or the IFOF. Since the influential work of Andreas Vesalius (De humani corporis fabrica, 1543, Padua Italy) anatomists are drawing, dissecting and reconstructing WM pathways with increasing precision as technology advances. It would not be unexpected if the SAF could reemerge from textbooks of human brain anatomy in a similar way.

#### CONCLUSION

We have proposed the existence of the SAF in the human brain, a new WM bundle identified with dMRI based FT. We showed that it is consistent within several cohorts, across different acquisition settings and processing strategies. However, there are still uncertainties about the true underlying anatomy of the structure as evidenced by our complementary PLI and brain dissection findings.

### DATA AVAILABILITY

All datasets generated for this study are included in the manuscript and the supplementary files.

### ETHICS STATEMENT

#### Dissection (Paris)

The study was approved by the local ethics committee of Pôle Neurosciences of Lariboisière Hospital. This study complied with the ethical standards.

#### REFERENCES

Alexander, D. C., Dyrby, T. B., Nilsson, M., and Zhang, H. (2017). Imaging brain microstructure with diffusion MRI: practicality and applications. NMR Biomed. doi: 10.1002/nbm.3841 [Epub ahead of print].

### Dissection (Trento)

Dissections were performed in the context of the Structural and Functional Connectivity Lab (SFC-Lab) Project, regularly approved by the Ethical Commitee of the APSS di Trento (authorization N◦ 1837 released on September 26th 2013 and renewed with amendment on September 7th 2018).

### PLI

The brains were taken from persons without history of neurologic or psychiatric disease, who donated their body for anatomical study. All brains were collected from the body donor program of the Institute of Anatomy at the Technical University Aachen (RWTH).

### AUTHOR CONTRIBUTIONS

SD, AH, and AL: protocol design and data analysis. SD, AH, FCorr, MTS, SS, FCors, AB, LP, EM, and AL: data interpretation. FCorr, MTS, SS, FCors, AB, LP, and EM: neuroanatomy. HA: PLI. SD, AH, FCorr, MTS, SS, FCors, AB, LP, EM, HA, TP, and AL: manuscript preparation. AH, FCorr, MTS, SS, FCors, AB, LP, MV, DJ, EM, HA, JE, TP, and AL: manuscript revision. AL: principal investigator.

### FUNDING

The research of AL, AH and SD is supported by VIDI Grant 639.072.411 from the Netherlands Organization for Scientific Research (NWO). The UK Medical Research Council and the Wellcome Trust (Grant ref: 092731) and the University of Bristol provided core support for ALSPAC. Acquisition of MRI data was funded by the National Institutes of Health (R01MH085772 to TP).

### ACKNOWLEDGMENTS

We are extremely grateful to the families who took part in this study and the ALSPAC team, including midwifes, interviewers, computer and laboratory technicians, clerical workers, research scientists, volunteers, managers, receptionists and nurses. A pre-print of this work has been released at bioRxiv (David et al., 2018).

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnana. 2019.00024/full#supplementary-material

Andersson, J. L. R., Graham, M. S., Zsoldos, E., and Sotiropoulos, S. N. (2016). Incorporating outlier detection and replacement into a non-parametric framework for movement and distortion correction of diffusion MR

Anderson, A. W. (2005). Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magn. Reson. Med. 54, 1194–1206. doi: 10.1002/mrm.20667

images. Neuroimage 141, 556–572. doi: 10.1016/j.neuroimage.2016. 06.058


of diffusion magnetic resonance parameters. Magn. Reson. Med. 73, 2174–2184. doi: 10.1002/mrm.25351


**Conflict of Interest Statement**: 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.

Copyright © 2019 David, Heemskerk, Corrivetti, Thiebaut de Schotten, Sarubbo, Corsini, De Benedictis, Petit, Viergever, Jones, Mandonnet, Axer, Evans, Paus and Leemans. 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.

# A Missing Connection: A Review of the Macrostructural Anatomy and Tractography of the Acoustic Radiation

#### Chiara Maffei1,2 \*, Silvio Sarubbo<sup>3</sup> and Jorge Jovicich2,4

<sup>1</sup> Harvard Medical School, Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Charlestown, MA, United States, <sup>2</sup> Center for Mind/Brain Sciences - CIMeC, University of Trento, Trento, Italy, <sup>3</sup> Division of Neurosurgery, Structural and Functional Connectivity Lab Project, S. Chiara Hospital, Trento Azienda Provinciale per i Servizi Sanitari (APSS), Trento, Italy, <sup>4</sup> Department of Psychology and Cognitive Sciences, University of Trento, Trento, Italy

The auditory system of mammals is dedicated to encoding, elaborating and transporting acoustic information from the auditory nerve to the auditory cortex. The acoustic radiation (AR) constitutes the thalamo-cortical projection of this system, conveying the auditory signals from the medial geniculate nucleus (MGN) of the thalamus to the transverse temporal gyrus on the superior temporal lobe. While representing one of the major sensory pathways of the primate brain, the currently available anatomical information of this white matter bundle is quite limited in humans, thus constituting a notable omission in clinical and general studies on auditory processing and language perception. Tracing procedures in humans have restricted applications, and the in vivo reconstruction of this bundle using diffusion tractography techniques remains challenging. Hence, a more accurate and reliable reconstruction of the AR is necessary for understanding the neurobiological substrates supporting audition and language processing mechanisms in both health and disease. This review aims to unite available information on the macroscopic anatomy and topography of the AR in humans and non-human primates. Particular attention is brought to the anatomical characteristics that make this bundle difficult to reconstruct using non-invasive techniques, such as diffusion-based tractography. Open questions in the field and possible future research directions are discussed.

Keywords: acoustic radiation, auditory system, sensory pathways, auditory pathways, auditory tract, diffusionbased tractography

**Abbreviations:** AC, auditory cortex; AR, acoustic radiation; CC, corpus callosum; CR, corona radiata; CSD, constrained spherical deconvolution; dMRI, diffusion magnetic resonance imaging; DWI, diffusion weighted imaging; EC, external capsule; EEG, electroencephalography; FA, fractional anisotropy; FS, sylvian fissure; HC, healthy cotrols; HG, Heschl's gyrus; Ic, inferior culliculus; IC, internal capsule; ICPL, posterior limb of the internal capsule; ICSL, sub-lenticular part of the internal capsule; ILF, inferior longitudinal fasciculus; LGN, lateral geniculate nucleus; MEG, magnetoencephalography; MGN, medial geniculate nucleus; ODF, orientation distribution function; OR, optic radiation; ROI, region of interest; PAC, primary auditory cortex; SNR, signal-to-noise ratio; STG, superior temporal gyrus; STP, superior temporal plane; TH, thalamus; WM, white matter; XRT, X-ray therapy.

#### Edited by:

Ricardo Insausti, University of Castilla La Mancha, Spain

#### Reviewed by:

David Reser, Monash University, Australia Hisayuki Ojima, Tokyo Medical and Dental University, Japan

#### \*Correspondence:

Chiara Maffei cmaffei@mgh.harvard.edu

Received: 04 April 2018 Accepted: 15 February 2019 Published: 07 March 2019

#### Citation:

Maffei C, Sarubbo S and Jovicich J (2019) A Missing Connection: A Review of the Macrostructural Anatomy and Tractography of the Acoustic Radiation. Front. Neuroanat. 13:27. doi: 10.3389/fnana.2019.00027

### INTRODUCTION

fnana-13-00027 March 5, 2019 Time: 19:7 # 2

The acoustic radiation (AR) represents a highly-myelinated group of axonal projections and constitutes one of the primary sensory pathways of the primate brain, carrying auditory information from the thalamus to the cortex. The connectivity pattern of these fibers has been described in some detail in cytoarchitectonic and myeloarchitectonic studies of non-human primates (Polyak, 1932; Mesulam and Pandya, 1973; Morel et al., 1993; Hackett et al., 1998) and, at a more macroscopic level, in a few histological studies in humans (Flechsig, 1920; Pfeifer, 1920; Rademacher et al., 2002; Bürgel et al., 2006). However, the information obtained from non-human primate studies cannot be transferred directly to the human brain. Furthermore, such studies have focused mostly on the cytoarchitectonic aspects of the auditory cortices and their intrinsic connectivity, with little emphasis on the anatomical course of the AR itself. In humans, tracing studies are impossible in vivo and have restricted applications in post-mortem brains (Mesulam, 1979; Tardif and Clarke, 2001), while limited information can be drawn from old myeloarchitectonical post-mortem studies.

The advent of diffusion magnetic resonance imaging (dMRI) (Basser et al., 1994) and tractography (Mori et al., 1999) has made it possible to investigate the anatomy of the major white matter (WM) bundles of the human brain in vivo and non-invasively (Catani et al., 2002; Catani and de Schotten, 2008; Lawes et al., 2008). However, the AR constitutes a notable exception in this sense. This primary sensory bundle is largely absent from most tractography studies investigating audition and language and from human WM atlases (Thiebaut de Schotten et al., 2011). This is mainly due to the intrinsic anatomical characteristics of these fibers, which go beyond the current limits of dMRI tractography methods (Behrens et al., 2007; Jones and Cercignani, 2010; Daducci et al., 2016). Therefore, the diffusion-based tractography reconstruction of the AR remains highly challenging at present, discouraging its in vivo anatomical investigation in humans.

However, overcoming or circumventing these methodological considerations is essential, as successful in vivo reconstruction of the human auditory tract is of great importance for both clinical applications (e.g., pre-surgical mapping) and basic neurobiological research. Reliably revealing the 3D characteristics of this tract would help in correlating the anatomical and functional aspects of audition and in the study of human-only cognitive functions, such as language, both in healthy and pathological conditions.

The main aim of this review is to emphasize the paramount need for characterizing the human AR in both clinical and scientific contexts, which has yet to be done for a number of reasons that we discuss. This review collates the available information from primate studies on the anatomy, topography, and course of the AR, with particular emphasis on the anatomical features that make this tract extremely challenging to study, even for state-of-the-art dMRI tractography techniques. Additionally, recent attempts to reconstruct the AR using diffusion-based tractography methods will be discussed. Finally, open questions in the field will be presented and possible future research directions considered.

### THE ACOUSTIC RADIATION IN PRIMATES

The auditory system of mammals is a complex network of parallel and overlapping axonal projections that connect subcortical nuclei and cortical regions. It encodes and transmits stimuli coming from the acoustic environment, enabling an organism to detect a sound in its environment, determine the direction from which it originated, discriminate among potential sources, and thereby, react or communicate with conspecifics. **Figure 1** provides a schematic representation of the human ascending auditory system and its sub-cortical relays, as it is commonly described in both the scientific literature and neuroanatomy textbooks (Moller, 2006; Brugge, 2013). The thalamo-cortical projections of this system, which connects the medial geniculate nucleus (MGN) to auditory cortex, constitute the AR (**Figure 1**).

Anatomical knowledge of the human AR mainly comes from pioneering investigations at the beginning of the 20th century (Dejerine and Dejerine-Klumpke, 1895; Flechsig, 1920; Pfeifer, 1920). Since these early studies, very limited additional information on this structure has been reported for humans (Rademacher et al., 2001, 2002; Bürgel et al., 2006). Most of the anatomical and functional organization of the mammalian auditory system has been inferred from animal studies, mainly non-human primates and cats (Morel and Kaas, 1992; Hashikawa et al., 1995; Hackett et al., 1998; Kaas and Hackett, 2000; Hackett et al., 2001; Jones, 2003; de la Mothe et al., 2006; Lee and Winer, 2008). However, these studies primarily focus on the topographical mapping between the thalamus and auditory cortex without documenting the spatial connection pattern and course of these fibers (Jones, 2003; Hackett, 2011). As a consequence, most neuroanatomical books only report schematic drawings of this pathway (e.g., **Figure 1**). Furthermore, although comparative studies have shown similar features across human and non-human primate brains, cortical and subcortical architectonic differences, as well as cognitive dissimilarities,

extra-lemniscal pathway by dashed lines. CN, cochlear nucleus; SOC, superior olivary complex; LL, lateral lemniscus; Ic, inferior culliculus; MGN, medial geniculate nucleus; AR, acoustic radiation.

exist (Galaburda and Sanides, 1980; Galaburda and Pandya, 1983; Schmahmann et al., 2007; Passingham, 2009; Thiebaut de Schotten et al., 2012). Thus, it is important to understand the commonalities among primate species and to identify potentially unique aspects of the human auditory system that might be related to our ability to perceive and process languagespecific stimuli. In vivo diffusion imaging techniques constitute a powerful tool in investigating these topics.

In the next section, we briefly review the AR microstructural topography, as described in animal studies, and its macrostructural anatomy, which has been gleaned from human research. Our goal is to provide a more complete anatomical profile of this bundle across species and to highlight the existing gap of topographical information between invasive and non-invasive studies. In particular, we focus on the AR in vivo imaging literature, reviewing recent attempts to visualize this tract using tractography techniques.

### The Acoustic Radiation in Invasive Studies

The axonal connections between the thalamus and the auditory cortex have been investigated in animals using different invasive techniques. These fibers stem from the medial geniculate nucleus of the thalamus (MGN), as first described by von Monakow (1882), and contact a specific area on the posterior part of the Sylvian fissure (Minkowski, 1923; Polyak, 1932), which has been described as a "rudimentary transverse temporal gyrus" (Walker, 1937) in monkeys. In most species, three major divisions of the MGN are identified: ventral (or principal), dorsal (or posterior), and medial (or magnocellular) (Winer et al., 2001; Jones, 2003). Each of these divisions has unique connections to different cortical regions. The ventral division receives input from the central nucleus of the inferior culliculus (Ic) and almost exclusively projects to what is defined as the core region of the auditory cortex (Mesulam and Pandya, 1973; Burton and Jones, 1976; Morel and Kaas, 1992; Morel et al., 1993; Hashikawa et al., 1995; Rauschecker et al., 1997). This region is distinguished by dense immunoreactivity for the calciumbinding protein parvalbumin, as most of its inputs come from the ventral MGN parvalbumin immunoreactive cells (Molinari et al., 1995), even if some connections with the other MGN divisions appear to exist (Luethke et al., 1989; Morel et al., 1993). This region occupies a portion of the caudal superior temporal plane (postero-medial part of the Heschl's gyrus in humans) and it is characterized by a dense population of small granule cells (e.g., koniocortex) with a well-developed layer IV (Merzenich and Brugge, 1973; Seldon, 1981). Both the ventral MGN and the core region show a tonotopical organization in which the representation of frequencies is spatially organized. This suggests a topographical organization of fibers connecting similar frequency domains in these two structures (Burton and Jones, 1976; Molinari et al., 1995). The core region is surrounded by a secondary narrow belt region and a third, more lateral region that occupies the lateral surface of the superior temporal gyrus (Hackett et al., 1998; Kaas and Hackett, 2000). This latter "parabelt" region is generally considered to be a higher-order auditory region or auditory association cortex that integrates auditory with non-auditory multisensory information (Hackett et al., 1998). These regions are less responsive to pure tone sounds, preferring more complex sounds, and do not show the clear tonotopical organization typical of the core region (Rauschecker et al., 1995; Jones, 2003). The dorsal and medial divisions of the MGN constitute the major inputs to these secondary auditory association regions. These nuclei receive inputs from the external nucleus of the Ic, as well as from lower brainstem relays, and bypass the core region to project to secondary auditory and other cortical regions (Rauschecker et al., 1997; Jones, 2003; Winer and Lee, 2007).

Projections from the ventral MGN to the core region correspond to the most direct classic auditory pathway, also called the lemniscal pathway. These parallel connections are tonotopically organized and their neurons show sharp responses to tones (Morel et al., 1993). Direct projections from the other MGN divisions to secondary auditory cortical regions are part of the extralemniscal (or non-lemniscal) auditory pathway. These fibers are separate from but lie adjacent to those of the lemniscal pathway in the ascending auditory system. In addition, their subdivision continues at the cortical level, where the non-lemniscal pathway has stronger and more diffuse connections with regions surrounding the core region (Lee and Sherman, 2010). These pathways are less tonotopically organized and their neurons demonstrate fewer sharp responses to sounds.

Classical studies in non-human primates focus on the topography of the connectivity between the MGN and the auditory projection cortical territory but provide no information on either the course and extension of the AR tract itself or its relationship with the other WM pathways of the brain. Using Marchi axonal degeneration, Polyak (1932) provided a detailed description of the course of this tract in rhesus macaques (**Figure 2**). He describes a dense bundle of closely assembled fibers that leaves the MGN, turns laterally and crosses the most ventral portion of the internal capsule (IC) immediately above the lateral geniculate nucleus (LGN) of the thalamus (**Figure 2**). At this level, these fibers are distinguishable from somatosensory fibers because of their nearly horizontal orientation. The AR then bends ventrally and reaches the external capsule (EC) by passing through the ventral edge of the posterior putamen. Once there it meets other projection and association bundles before finally reaching the WM of the superior temporal convolution close to the Sylvian fissure (**Figure 2**). He describes the AR as a regularly arranged projection system, where fibers lie parallel to one another until gradually diverging only when approaching the cortex.

Classical topographical descriptions in humans (**Figure 3**) provide a very similar description (Dejerine and Dejerine-Klumpke, 1895; Flechsig, 1920; Pfeifer, 1920) of the AR, which may be summarized as follows. The AR leaves the MGN and travels in an antero-lateral direction; it then passes through the posterior portion of the IC, proceeds along the corona radiata and curves around the inferior portion of the circular sulcus of the insula before entering the transverse temporal gyrus of Heschl (HG) in a ventral-to-dorsal direction (Pfeifer, 1920).

At this macrostructural level, both animal and human studies delineate a bundle with a transverse orientation that lies adjacent to and crosses over other main WM bundles before reaching the auditory cortex. In both animals and humans, the AR intermingles with the fibers of the IC in its most posterior portion. Furthermore, Dejerine and Dejerine-Klumpke (1895) described the close proximity of the AR to the optic radiation (OR) at the stemming point in the thalamus, defining this region as the "carrefour sensitive" (sensory intersection). Polyak (1932) also studied this region, stating that the OR and AR, although lying close together, are completely separate: the AR is located more anterior and crosses at a right angle above the OR, which follows a posterior direction in the sagittal plane. As will be discussed in the following section, this configuration and certain other anatomical features of the AR pose serious challenges to its 3D tractography reconstruction.

The myeloarchitectonic maps from Dejerine and Dejerine-Klumpke (1895) and Flechsig (1920), while being of invaluable historical significance, cannot be used to extract precise anatomical information that can be applied to modern brain atlases or neuroimaging studies. More recently, radiological information about the AR anatomical organization was obtained in human post mortem myelin-stained sections (Rademacher et al., 2002; Bürgel et al., 2006). These studies confirm the classical topographical description of the acoustic fibers and are of great importance as they represent the main reference framework for in vivo imaging studies of this brain region and provide the opportunity to investigate inter-subject variability and hemispheric asymmetry. Previous studies have found that the AR does not enter the lenticular nucleus, but rather runs

human brain (coronal view). (B) Coronal cut through the middle section of the thalamus. (C) Axial section of the brain; cut through the inferior thalamus. In the three panels, fibers belonging to the acoustic radiation have been highlighted in green and fibers of the internal capsule in pink. The acoustic radiation projects from the thalamus (TH) to the first temporal circonvolution (T1) and passes through the sub-lenticular and posterior segment of the internal capsule. This map clearly highlights the crossing between these two fiber systems. CC, Corpus callosum; EC, external capsule; ICPL, Posterior limb of the internal capsule; ICSL, sub-lenticular part of the internal capsule; CR, corona radiata; LGN, lateral geniculate nucleus; MGN, medial geniculate nucleus; T1, first temporal circonvolution; TH, thalamus (adapted from Dejerine and Dejerine-Klumpke, 1895; https://archive.org, public domain).

dorsally to the OR, crossing the temporal isthmus as it ascends to the auditory cortex (Pfeifer, 1920; Polyak, 1932; Bürgel et al., 2006). Fanning of these fibers in HG creates a hat-like structure that covers the posterior end of the lenticular nucleus (Rademacher et al., 2002).

According to Flechsig (1920), these fibers are divisible into two bundles, one of which ascends near the Ec and enters the auditory cortex from the superior-posterior side. The other bundle courses for some distance in the company of the OR before passing behind and below the fossa sylvii where it pierces the bases of the middle and inferior temporal gyri to reach the transverse temporal gyrus or gyri. Similarly, early evidence from animal studies suggests that the AR is subdivided into dorsal and ventral components (von Monakow, 1882), although, Rademacher et al. (2002) found more recently only a single and heavily myelinated bundle.

Overall, the macrostructural description of the AR in humans resembles the description of this tract in non-human primates. The origin and termination of this bundle, together with the extension and relationship to other WM bundles, is maintained. However, the macro- and micro-anatomical correspondence between the different cortical regions in monkeys and humans is not straightforward (Baumann et al., 2013). Compared to nonhuman primates, the human cortical surface of the auditory regions demonstrates additional gyri and higher inter-subject and interhemispheric variability (Galaburda et al., 1978; Hackett et al., 2001), both of which may affect the AR anatomy. The human auditory cortex also shows higher differentiation into

sub-regions as compared to non-human primates, and the core region is larger than the belt region (the opposite is true for monkeys) (Fullerton and Pandya, 2007). Both the differences in cortical anatomy between non-human primates and humans and the major variability of cortical and subcortical structures across subjects and hemispheres in humans (Bürgel et al., 2006) raise interesting questions about the possibility of a relationship between such morphological differences and human-only language abilities.

### The Acoustic Radiation in Non-invasive Tractography Studies

Diffusion MRI (dMRI) tractography allows for the investigation of WM architecture in the human brain non-invasively in vivo. Since its first applications (Mori et al., 1999), most of the well-known WM bundles of the human brain have been reconstructed using diffusion-based tractography methods (Catani and Thiebaut de Schotten, 2008; Lawes et al., 2008). Despite its potentials, dMRI tractography has several important limitations that have been discussed in depth in the literature (e.g., Jones and Cercignani, 2010; Thomas et al., 2014; Maier-Hein et al., 2017). In principle, these limitations affect most tractography reconstructions, but here we focus on how they particularly affect the 3D reconstruction of the AR. Reconstructing the AR three-dimensionally is highly challenging at present due to the anatomical features described in the previous sections: its relatively small size, transversal orientation, and location in a region with a high density of crossing fibers. This largely prevents the inclusion of this particular tract in most tractography investigations.

As shown in the previous section (see **Figures 2**, **3**), in its medio-lateral course from the MGN to the HG, the AR lies in a nearly horizontal position and, for this reason, crosses some of the major fiber systems of the human brain: internal capsule, external capsule, and posterior thalamic radiation (Maffei et al., 2018). Resolving the fiber crossing is a wellknown challenge in dMRI (Tuch et al., 2002; Dell'Acqua and Catani, 2012; Jeurissen et al., 2013). The classic tensor model (Basser et al., 1994) is capable of characterizing only one main fiber orientation per voxel and it has been shown to constantly fail in regions where voxels contain complex fiber architectures (Behrens et al., 2007; Jbabdi and Johansen-Berg, 2011). The impact of this limitation is particularly evident for non-dominant tracts, given that the orientation produced by the tensor will be closest to the largest contributing direction in most cases. This effect is amplified at the low resolution of commonly available diffusion protocols due to within-voxel partial volume averaging effects (Tournier et al., 2011). When implemented in tractography studies, the diffusion tensor model has proven to be incapable of detecting the 3D profile of the AR. Streamlines are either truncated when entering voxels containing major inferior-superior orientations or erroneously embedded in the reconstruction of these major projection bundles, with no visible streamlines contacting the HG (Behrens et al., 2007; Crippa et al., 2010; Berman et al., 2013). This has likely been the primary factor preventing the investigation of the auditory system by means of diffusion-based tractography. Some studies have used the diffusion tensor to investigate the WM microstructure of the auditory system, limiting the structural investigation of the auditory pathways to the extraction of mean quantitative diffusion measures [e.g., fractional anisotropy (FA)] from specific regions of interest (ROI) (Chang et al., 2004; Lee et al., 2007; Lin et al., 2008; Wu et al., 2009). However ROIbased analysis can lead to inaccurate results, especially for WM tracts that are extremely variable across subjects, such as the AR (Rademacher et al., 2002). Therefore, it is typically preferable to map the exact anatomy of such WM tracts in individual subjects/patients.

To address the intrinsic limitations of the tensor formalization, more advanced models have been introduced that can better account for fibers crossing, by modeling more than one fiber population per voxel (Tuch, 2004; Tournier et al., 2008; Descoteaux et al., 2009). These models open up the possibility of propagating streamlines through crossing fiber regions, thus allowing the reconstruction of non-dominant WM bundles, such as the AR. However, together with the low-level diffusion model employed, other parameters play a role in the accurate reconstruction of WM bundles, such as the tractography parameters chosen and the strategy to define inclusion ROIs. Here we report studies that reconstruct the AR 3D tractography profile in vivo using multi-fiber-based models (**Figure 4**), taking into consideration the different combinations of acquisition parameters, diffusion models and ROIs selection strategies that they used (**Table 1**).

Behrens et al. (2007) were able to visualize the course of the AR from the MGN to the cortex using the ball-and-stick model and probabilistic tractography (Behrens et al., 2003) (**Figure 4A**). The authors demonstrate how this multi-fiber model can reconstruct the AR, overcoming the limitations of the tensor model. Three following studies were then able to demonstrate the reliability of this model for successful in vivo reconstruction of the profile of the auditory tract in both healthy subjects and those with tinnitus (Crippa et al., 2010; Javad et al., 2014; Profant et al., 2014). In these studies, different combinations of inclusion ROIs were used to isolate the AR, including the Ic, the MGN, and both functionally and manually defined HG. Nevertheless, the profile of the 3D tractography reconstruction looks visually similar across the four studies and shows connections between the posterior thalamus and the auditory region WM. However, artifacts are visible along the inferior-superior axis in the middle part of the AR at the level of the crossing with the IC and these false positive reconstructions are likely to be related to the probabilistic nature of the diffusion model and tractography algorithm used. Despite the inclusion of false positive signals in the reconstructions, the ball-and-stick model has the important advantage of being accessible to low b-value (b = 1000 s/mm<sup>2</sup> ) diffusion protocols which makes it suitable for clinical investigations of the AR. Berman et al. (2013) (**Figure 4B**) used a solid-angle q-ball model (Aganj et al., 2010) and probabilistic tractography to successfully reconstruct the auditory connections between the MGN and HG. This reconstruction looks anatomically very accurate and free of false positive artifacts. However, q-ball methods require higher b-value diffusion data (b ≥ 3000 s/mm<sup>2</sup> ). In these models,

FIGURE 4 | The tractography reconstruction of the acoustic radiation (AR) in three different studies. (A) 3D volumetric reconstruction of the right AR in one subject using a multi-tensor model and probabilistic tractography. Voxels are color-coded from 10 (red) to 50 (yellow) samples passing through the voxel (adapted with permission from Behrens et al., 2007). (B) The panel shows the streamlines of the right AR in axial view, as reconstructed in one subject using q-ball imaging and probabilistic tractography. The location of the thalamus and auditory cortex (AC) are specified by the blue boxes. The figure also shows the orientation distribution functions (ODF) corresponding to the inferior longitudinal fasciculus (ILF), highlighted by the green arrow (adapted with permissions from Berman et al., 2013). (C) Left: 3D tractography reconstruction of the right AR in one subject using constrained spherical deconvolution models and probabilistic tractography. A 3D rendering of the thalamus is also shown in gray and the borders of HG in green. Right: Klinger's post-mortem blunt dissection of the right AR (modified from Maffei et al., 2018). In all three panels the location of thalamus (blue arrow) and auditory cortex (white arrow) are highlighted.

the angular resolution of the reconstructed diffusion profiles is increased and the crossing fiber configurations are correctly represented (Tournier et al., 2013), although they pose limitations to its use in clinical populations. Our group (Maffei et al., 2018) used ultra-high b-value Human Connectome Project diffusion data-sets (Fan et al., 2015) and spherical deconvolution (Tournier et al., 2008) to reconstruct AR streamlines using probabilistic tractography. We compared the results to Klinger's post-mortem blunt micro-dissections (**Figure 4C**), a method based on a brain freezing technique optimized to reveal WM (Ludwig and Klingler, 1956). This approach has been used in several studies to evaluate tractography accuracy (Fernaìndez-Miranda et al., 2015; De Benedictis et al., 2016; Pascalau et al., 2018). The obtained reconstructions agreed with the AR anatomy revealed in the post-mortem dissections and no additional exclusion regions were needed to isolate the AR profile. However, the ultra-high b-values used in this study (b ≤ 10,000 s/mm<sup>2</sup> ) are very rarely achievable, even in research settings. More recently, Yeh et al. (2018) published a population-averaged atlas of several WM connections, including the AR, using q-space diffeomorphic reconstruction (QSDR) (Yeh and Tseng, 2011) and deterministic tractography for multiple fiber orientations. In their approach, the high angular and spatial resolution of the data (1.25 mm isotropic, and b ≤ 3000 s/mm<sup>2</sup> ) and the large sample (842 subjects) allowed them to reduce the rate of false positive artifacts. However, the AR profile shown in this study, while correctly originating at the posterior thalamus, does not reach the expected auditory cortex on the superior side of the temporal lobe.

Overall, the reconstructed 3D profiles shown in these studies are in accordance with the macrostructural landmarks defined by classic anatomical studies: streamlines originate in the posterior thalamus and course in an antero-lateral direction to terminate in the temporal lobe (Dejerine and Dejerine-Klumpke, 1895). However, AR reconstructions are still highly variable across studies. In particular, while showing similar profiles at the thalamic stemming region, reconstructions differ as they approach the cortex, either falling short of reaching the HG (Crippa et al., 2010; Yeh et al., 2018) or creating false positive artifacts at the intersection with vertically oriented fibers (Javad et al., 2014; Maffei et al., 2018). Moreover, disagreement about the relationship with neighboring tracts exists. Berman et al. (2013) suggest that AR streamlines cross the inferior longitudinal fasciculus (ILF), while Behrens et al. (2007) and Javad et al. (2014) claim that they cross the OR. In a recent work from our group, we did not find that the AR is in close proximity to the ILF (Maffei et al., 2018), supporting older studies that report no crossing between the AR and OR (Pfeifer, 1920; Polyak, 1932), as described above. In addition to variability across studies, low reproducibility across subjects is also reported. Some groups have been able to reconstruct AR tracts on both hemispheres on 100% of subjects (Behrens et al., 2007; Profant et al., 2014; Maffei et al., 2017). In contrast, even when using similar diffusion models (e.g., ball-and-stick), other studies report reconstructions successful in both hemispheres in much lower proportions, such as 35–50% (Crippa et al., 2010) or 71–86% (Javad et al., 2014).

We suggest that the present variability and low reproducibility in the reconstructed AR profile are related to a combination of some of the specific characteristics of the AR that make its tractographic reconstruction quite challenging, even for state of the art tractography techniques: its anatomical location, small size, and inter-individual anatomical variability.

As alluded to in the previous section, the AR constitutes a compact but relatively short and small bundle that lies horizontally in a region with a high density of vertical fibers. Even if multi-fiber models proved capable of representing this crossing, the degree to which this crossing can be accurately resolved in the final 3D reconstruction also depends on the tractography algorithm used and the intrinsic angular resolution of the dMRI data (Tournier et al., 2013). For example, using a higher b-value (Berman et al., 2013; Maffei et al., 2018) might help improve accuracy of results relative to those obtained with lower b-value data (Crippa et al., 2010). However, in these studies several exclusion ROI have been employed to either constrain tractography (Behrens et al., 2007) or clean the results (Crippa et al., 2010). This renders the accuracy of


TABLE 1 | The table reports the main acquisition and tractography parameters used to reconstruct the AR in the listed studies.

The acquisition parameters column provides the number of diffusion encoding directions, the b-factor (s/mm<sup>2</sup> ), and the spatial resolution (mm). The inclusion ROIs column indicates whether seed and target regions were selected manually or with other methods. Some tractography studies performed reconstructions in both directions and results were summed. The success rate column reports the percentage of subjects showing successful AR tractography reconstructions in both hemispheres. DWI, diffusion MRI encoding directions; CG, Cingulate Gyrus; CSD, Constrained Spherical Deconvolution; HG, Heschl's Gyrus (gray and white matter); HG-WM, Heschl's Gyrus (only white matter); Ic, Inferior Culliculus; MGN, Medial Geniculate Nucleus; QSDR, q-space diffeomorphic reconstruction; TH, Thalamus; HC, healthy controls.

the final reconstructions sensitive to the selection of these ROI, complicating comparisons across studies.

The use of different ROI selection strategies can strongly affect the resulting tractography reconstructions. The HG is a complex structure that shows large variability in sulcal landmarks across subjects and hemispheres (Rademacher et al., 2001), whereas the MGN is a very small structure, varying from 74 to 183 mm<sup>3</sup> (Kitajima et al., 2015), making it difficult to locate in neuroimaging data and also highly variable across individuals and hemispheres (Rademacher et al., 2002). Therefore, this variability poses difficulties in the selection of the ROIs used to initiate the tractography reconstruction or to perform the virtual dissections. For example, while reliable automatic segmentation tools for the entire thalamus are available in different public software packages (e.g., FSL<sup>1</sup> , Freesurfer<sup>2</sup> ), it is far more challenging to automatically segment smaller structures, such as the MGN, due to both their size and their lower MRI contrast with neighboring WM. Alternative solutions exist but are also challenging. For example, subject-specific manual segmentations of MGN can lead to high anatomical accuracy, but they are very time intensive and, therefore, costly to do in large subject cohorts. Brain atlases also can be used to define the MGN, however given the small size and inter-subject variability, atlas-driven segmentations are unlikely to provide a good anatomical match for all subjects. The development of more accurate automatic parcellation techniques for the thalamic nuclei is expected to improve the accuracy of seed-to-target definition and, thus, of the resulting tractography reconstructions. Recently, a new and promising probabilistic atlas of the thalamic nuclei has been proposed based on a combination of ex vivo MRI and histology (Iglesias et al., 2018).

Future research in the field should focus on how to improve the sensitivity and reproducibility of the tractography reconstruction of this bundle. This could be achieved by inputting prior anatomical knowledge in the tractography process, as it is implemented in global-tractography-based frameworks (Yendiki et al., 2011) or more recently developed bundlespecific algorithms (Rheault et al., 2019). Parallel to these advances at the tractography level, there have been efforts in validating tractography results at a micro-anatomical scale. As this validation process progresses, we expect to expand our knowledge of the exact boundaries of the human AR and, consequently, better inform its tractography reconstruction and improve accuracy of tractography results. A more accurate tractography investigation of the AR could expand our structural and functional knowledge of the auditory system, as proposed in the next section.

<sup>1</sup>https://fsl.fmrib.ox.ac.uk/fsl/

<sup>2</sup>http://surfer.nmr.mgh.harvard.edu/

The reliable in vivo reconstruction of the AR in humans may help the exploration of the neuro-anatomical and functional mechanisms underlying auditory processing and language comprehension. The precise characterization of the AR can provide information useful for clinical applications, such as in diagnosis and treatment of hearing and speech disorders, recovery from injury, and performance of interventions that can damage the AR, such as brain surgery or radiation treatments. This section provides a brief review of basic and clinical research areas that could benefit from an improved characterization of the AR.

#### Language and Auditory Perception

The ability to communicate through speech is quintessentially human. However, the anatomical organization and the functional mechanisms underlying speech comprehension in the brain are still not understood completely. The acoustic information that reaches the primary auditory cortex via the AR fibers is processed within neural networks that depend on cortico-cortical shortand long-range connections involving temporal, parietal and frontal regions, as schematized in the dual-stream model (Hickok and Poeppel, 2007; Saur et al., 2008; Friederici, 2009). Within this processing network, it is unclear where language-specific processing starts and whether the auditory cortex is involved in speech-specific analysis. Some theories suggest that the left auditory cortex is specialized in processing temporal cues that are fundamental for speech comprehension (Zatorre et al., 2002; Poeppel, 2003) and that this language-specific encoding might actually start at the subcortical level (Hornickel et al., 2009). The diffusion-based reconstruction of the AR, and of the auditory pathways at large, could help address the structural-functional relationship of speech perception. At the structural level, it would be interesting to understand whether the AR exhibits a degree of leftward lateralization in its volume, as demonstrated for some of the other WM bundles implicated in language processing (Catani et al., 2007). Reports on the macroscopic volumetric asymmetry of cortical auditory regions have been known for some time (Von Economo and Horn, 1930; Galaburda et al., 1978; Geschwind and Galaburda, 1985; Penhune et al., 1996), but only one study specifically investigated the hemispheric lateralization of the AR (Bürgel et al., 2006). At the functional level, the recent association of tractography and neurophysiological techniques [such as magnetoencephalography (MEG) and electroencephalography (EEG)] opens interesting possibilities for investigating these topics. EEG metrics have been recently correlated with diffusion metrics in the investigation of the OR (Renauld et al., 2016). Similarly, EEG/MEG- and tractography-derived measures could be combined to investigate the relationship between temporal cortical regions and auditory function in both healthy subjects and patients.

On a finer scale, AR streamline terminations could be combined with functional MRI to provide critical insights into the subdivision of the auditory cortex, the borders of which are still not clearly defined using in vivo neuroimaging methods (Baumann et al., 2013). In this sense tractography could be used to investigate the topographical organization of the AR with respect to the different subdivisions of the auditory cortex. The different auditory cortical regions show a hierarchical organization in information processing (Tardif and Clarke, 2001), from highly specialized core regions to more integrated tertiary para-belt regions. This is confirmed by functional MRI studies, which suggest a gradient of increasingly more complex and abstracted processing from primary to higher-order auditory regions (Rauschecker et al., 1995; Humphries et al., 2014). Direct connections from MGN to secondary regions have been shown in monkeys (Rauschecker et al., 1997), suggesting that this functional organization might be maintained in the topographical organization of the AR fibers. At its present stage, tractography can locate and delineate the profile of major WM bundles with some accuracy, but it is very difficult to achieve precise site-to-site connectivity analysis with it; this limits the in vivo investigation of the WM topographical organization of the human brain. However, some studies use diffusion tractography to parcellate functionally different cortical regions (Rushworth et al., 2006; Anwander et al., 2007) and investigate the topographical organization of major bundles (Lee et al., 2016), and new methods have been proposed to advance the use of tractography for this purpose (Aydogan and Shi, 2016). In this scenario, it would be interesting to investigate whether some language-specific connections exist inside the AR and whether these project to higher-order languagespecific cortical regions. Moreover, this would also allow for the investigation of whether the tonotopical organization of the primary auditory cortex is reflected in its thalamocortical connections, as was recently shown in the mouse brain (Hackett et al., 2011). As dMRI acquisition (in particular, spatial resolution), diffusion modeling and tractography techniques improve, we will be able to bridge the gap between the micro-anatomical knowledge we have of the thalamocortical connections in animals and the macro-anatomical description in humans.

#### Language and Hearing Disorders

Damage to the auditory regions, most often the result of brain infarct or traumatic injury, has been associated generally with rare auditory syndromes, such as verbal auditory agnosia (Shivashankar et al., 2001), environmental auditory agnosia (Taniwaki et al., 2000), and cerebral (or central) deafness (Griffiths, 2002). However, until now, only one study investigated the extent of WM damage to the AR in a patient suffering from verbal auditory agnosia (Maffei et al., 2017). Investigating the extent of damage to the AR in patients suffering speech-related comprehension deficits would potentially enhance our understanding of the involvement of the AR in language processing.

Also, there is evidence that AR infarct can cause auditory hallucination (Woo et al., 2014), and that the extra-lemniscal pathway might be implicated in tinnitus perception (Moller et al., 1992). At present, studies investigating the auditory

pathways in these patients relied on WM ROI measurements (Lee et al., 2007; Lin et al., 2008), which only outline a portion of the underlying WM bundles and may not be representative of the entire tract. Being able to better understand the dynamics and location of such changes in the auditory pathways could help inform pathophysiological treatment strategies, such as repetitive transcranial stimulation (Langguth et al., 2010).

In congenitally and early deaf subjects, volumetric studies have outlined differences in gray and white matter of the auditory regions, compared to hearing subjects (Shibata, 2007; Kim et al., 2009), but to the best of our knowledge, no specific study on the AR in deaf subjects has been conducted to date. In addition to providing more detailed information on the anatomical changes occurring in the brain as a consequence of sensory deprivation, tractography of the AR may serve as an additional early diagnostic as well as complementary treatment tool in monitoring data in congenital hearing loss. This would help avoid delayed diagnosis that might lead to poor speech outcomes (Dedhia et al., 2018). Moreover, AR tractography reconstruction might be fundamental in assessing auditory pathway integrity before and after cochlear implantation, potentially predicting implant success (Huang et al., 2015).

The structural-functional relationship in language and hearing disorders can be further investigated by combining tractography reconstructions and more recently developed diffusion measures (Raffelt et al., 2012; Calamuneri et al., 2018) within both classical and more advanced tractography frameworks (Daducci et al., 2016). The application of these methods in clinical populations in the context of hearing disorders may help characterize axonal and myelination diseases (Ohno and Ikenaka, 2018) and auditory neuropathies (Moser and Starr, 2016; Ohno and Ikenaka, 2018).

Tractography reconstruction of the AR could help us investigate the anatomy of this tract in patients with hearing and/or language disorders, understand whether these fibers undergo structural reorganization in the case of auditory deprivation, and clarify the extent of AR damage in poststroke lesion profiles. This may be critical for shedding light on the functional-structural relationships of linguistic and nonlinguistic sound processing in the human brain.

#### Brain Surgical Planning

Investigating the functional and anatomical characteristics of the auditory fibers reaching the cortex, especially in relation to their implications for language function, would be important for surgical planning, such as in the case of tumor or epilepsy surgery (Wu et al., 2007; Farshidfar et al., 2014). The 3D reconstruction of major WM bundles is employed to plan and guide resections during surgery, and a functional atlas of human WM to drive well balanced onco-functional resections has been proposed recently (Sarubbo et al., 2015). In this context, diffusion-based virtual dissections have focused almost exclusively on language and sensory-motor structures (Chen et al., 2015). Possible reasons why AR fibers have not received much attention in the neurosurgical literature include the possibility that most of the non-linguistic auditory processing may happen at the brain-stem level and that auditory information is conveyed to both hemispheres, so that extensive bilateral damage is necessary for complete deafness (Griffiths, 2002). However, different sub-modal aspects of auditory processing, for example, those related to music perception and/or speech comprehension, might depend on the integrity of these projections (Hayashi and Hayashi, 2007; Baird et al., 2014). For cases of temporal lobe resection, the reliable virtual reconstruction of the AR might be critical for minimizing post-operative deficits in these domains. Also, it might serve in pre-operative assessments for cochlear implantation, as hearing recovery after implantation is influenced by the integrity of subcortical pathways (Vlastarakos et al., 2010).

#### Brain Radiation Oncology Planning

Today, X-ray therapy (XRT) is the standard of care for most brain tumors. However, XRT can damage normal brain tissue, causing neurocognitive deficits in different cognitive domains (Makale et al., 2016). The consideration of neuroimaging techniques for treatment planning is gaining importance as it helps avoid such complications by minimizing the unnecessary absorption of radiation in sensitive regions outside the tumor. Nevertheless, further improvement of radiation treatment will require tailored radiotherapy based on intra-treatment response (Wong et al., 2017).

Additionally, studies have shown XRT side effects in both gray (Bahrami et al., 2017) and white matter, although it is still largely unclear how variable sensitivity to radiation injury is across various regions of the brain (Connor et al., 2017). Kawasaki et al. (2017) show that tractography can help evaluate how radiation from XRT differentially affects WM regions and pathways. This knowledge, together with an understanding of how damage to such regions and pathways affects cognitive processes, could be used in the future to further optimize radiation treatment planning.

#### CONCLUSION

The anatomical and functional organization of the auditory system is still not well understood, particularly in humans. Successful in vivo tractographic reconstruction of the human auditory tracts is of great importance for clinical applications (e.g., pre-surgical mapping), as well as for basic research (e.g., language and auditory systems). This review outlines how the characterization of the AR has been limited by the methods used in the past and how advances in MRI acquisition and diffusion tractography methods offer the possibility to improve the characterization of this important WM tract. A few exciting potential research areas are suggested that would investigate anatomy and function concurrently in the same individual, both in health (e.g., the role of these tracts in language processing) and in disease (e.g., how the integrity of this tract relates to cognitive deficits). However, in order to obtain reliable reconstructions of the AR across subjects and protocols, additional work is needed to better understand how diffusion MRI acquisition

and tractography reconstruction strategies affect the AR 3D characterization and to validate tractography reconstructions at a micro-anatomical scale. Furthermore, as diffusion tractography is blind to the directionality of reconstructed fibers, the AR bundle could include both thalamo-cortical and corticothalamic projections, and therefore more studies are needed to differentiate between the afferent or efferent nature of these connections. Although these methodological challenges apply to diffusion MRI tractography in general, here we have focused on their relevance to the AR, a tract that has proven to be rather elusive for the reasons herein reviewed.

#### REFERENCES


### AUTHOR CONTRIBUTIONS

JJ conceived the review. SS revised the manuscript for neuroanatomical content. CM drafted the manuscript and designed the figures. All authors contributed to the final manuscript.

### ACKNOWLEDGMENTS

We thank Albert Galaburda for his insightful comments on a previous version of this manuscript.



and entorhinal pathways. Arch. Neurol. 36, 814–818. doi: 10.1001/archneur. 1979.00500490028004



neuronavigation. Neurosurgery 61, 935–949. doi: 10.1227/01.neu.0000303189. 80049.ab


**Conflict of Interest Statement:** 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.

Copyright © 2019 Maffei, Sarubbo and Jovicich. 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.

digital media

of impactful research

article's readership