An atlas of white matter anatomy, its variability, and reproducibility based on constrained spherical deconvolution of diffusion MRI

Virtual dissection of white matter (WM) using diffusion MRI tractography is confounded by its poor reproducibility. Despite the increased adoption of advanced reconstruction models, early region-of-interest driven protocols based on diffusion tensor imaging (DTI) remain the dominant reference for virtual dissection protocols. Here we bridge this gap by providing a comprehensive description of typical WM anatomy reconstructed using a reproducible automated subject-specific parcellation-based approach based on probabilistic constrained-spherical deconvolution (CSD) tractography. We complement this with a WM template in MNI space comprising 68 bundles, including all associated anatomical tract selection labels and associated automated workflows. Additionally, we demonstrate bundle inter- and intra-subject variability using 40 (20 test-retest) datasets from the human connectome project (HCP) and 5 sessions with varying b-values and number of b-shells from the single-subject Multiple Acquisitions for Standardization of Structural Imaging Validation and Evaluation (MASSIVE) dataset. The most reliably reconstructed bundles were the whole pyramidal tracts, primary corticospinal tracts, whole superior longitudinal fasciculi, frontal, parietal and occipital segments of the corpus callosum and middle cerebellar peduncles. More variability was found in less dense bundles, e.g., the fornix, dentato-rubro-thalamic tract (DRTT), and premotor pyramidal tract. Using the DRTT as an example, we show that this variability can be reduced by using a higher number of seeding attempts. Overall inter-session similarity was high for HCP test-retest data (median weighted-dice = 0.963, stdev = 0.201 and IQR = 0.099). Compared to the HCP-template bundles there was a high level of agreement for the HCP test-retest data (median weighted-dice = 0.747, stdev = 0.220 and IQR = 0.277) and for the MASSIVE data (median weighted-dice = 0.767, stdev = 0.255 and IQR = 0.338). In summary, this WM atlas provides an overview of the capabilities and limitations of automated subject-specific probabilistic CSD tractography for mapping white matter fasciculi in healthy adults. It will be most useful in applications requiring a reproducible parcellation-based dissection protocol, and as an educational resource for applied neuroimaging and clinical professionals.


Introduction
Characterizing the macroscopic structural organization of brain connectivity in vivo is central to understanding the human nervous system in health and disease. The advent of diffusion magnetic resonance imaging (dMRI) fiber tractography (FT) more than twenty years ago enabled significant progress in mapping major white matter (WM) fiber bundles described in anatomical and surgical literature. Initial work based on Diffusion Tensor Imaging (DTI) (Basser et al., 1994a(Basser et al., , 1994bMori et al., 1999) drove the development of white matter dissection protocols (Catani, 2006;Mori et al., 2009Mori et al., , 2008Wakana et al., 2007), which are still widely used today. This is because DTI is a simple and effective approach for reconstructing the core of large fasciculi (Catani et al., 2002;Stieltjes et al., 2001), data acquisition for DTI modelling requires a relatively short scan time (5-10 min) and standard pulse sequences are widely available owing to their regulatory approval for clinical practice. However, it is widely accepted that DTI suffers limitations that make it suboptimal for many tractography applications, particularly in a clinical setting (Farquharson et al., 2013;Mori and Tournier, 2014;Tournier et al., 2011). This has led to increased interest in more accurate and reliable approaches using high angular resolution imaging (HARDI) data (Bayrak et al., 2020;Bloy et al., 2012;Wasserthal et al., 2018;Yeh et al., 2018). At the same time, technological advances to accelerate data acquisition and reconstruction, such as multiband (Bouyagoub et al., 2020;Duan et al., 2015;Larkman et al., 2001, p. 1;Moeller et al., 2010) and compressed sensing (Lustig et al., 2008(Lustig et al., , 2007 in combination with improved computational efficiency, will increase the adoption of advanced reconstruction models in applications traditionally reserved for DTI. This means that there is a need for updated reference material based on more representative virtual dissections using HARDI data acquired on 3T scanners. Automated approaches relying on anatomical and orientational priors have been shown to considerably improve the accuracy of tract representations  and address the poor reproducibility, and operator dependency that confound the manual virtual dissection process (Kreilkamp et al., 2019;Maffei et al., 2021;Soares et al., 2013).There are currently only a limited number of approaches that automate the virtual dissection process from start to finish (Warrington et al., 2020;Wasserthal et al., 2018;Yendiki et al., 2011). However, often the underlying dissection protocol is not explicitly detailed or the workflow is based on model bundles defined a priori.
This work addresses both the currently unmet need for an updated HARDI human white matter atlas that is relevant for clinical research studies, and a standardized, reproducible virtual dissection approach based on anatomical definitions. Using CSD due to its versatility and potential application to clinical data (Calamuneri et al., 2018;Toselli et al., 2017;Wilkins et al., 2015) we extend earlier virtual dissection protocols based on DTI by providing a descriptive summary of the normal WM anatomy of 68 fiber bundles reconstructed using probabilistic tractography. We complement our anatomical descriptions with an open-source group atlas in MNI space and automated subject-specific virtual dissection software, which we call "Fun with Tracts" (FWT) that incorporates all the anatomical inclusion and exclusion labels per bundle. As probabilistic tractography produces inherently variable results, we also demonstrate how this may vary within and between individuals using open-source test-retest datasets from the human connectome project (HCP) (Van Essen et al., 2012) and the Multiple acquisitions for standardization of structural imaging validation and evaluation (MASSIVE) dataset . The result is a detailed and accessible reference for the virtual dissection of normal white matter anatomy using CSD tractography.

Summary
First, we describe and present an atlas of bundles generated with our literature-based dissection protocol (FWT) for 68 WM fasciculi based on the HCP-template data using a bundle-specific seeding and tractography approach described in more detail below. Tractogram images were generated using MRView (Tournier et al., 2019) and Surf-Ice (Rorden and Hanayik, 2022).
Second, we assessed inter-session/intra-subject variability by applying the FWT whole-brain tractography and bundle segmentation approach to the individual HCP test-retest datasets (Van Essen et al., 2012) and conducted similarity analyses between the resulting tractograms compared to the HCP-template bundles. Resulting weighted-dice similarity coefficient (wDSC) scores are shown in radar plots and were used to calculate single-rater agreement intra-class correlation (ICC) scores to highlight tractogram variability with respect to the template bundles. Additionally, we calculated descriptive statistics and generated violin plots for the wDSC scores resulting from pairwise comparisons of the output from each session (i.e., HCP scan 1 versus HCP scan 2). This provided a measure of FWT reproducibility within and between subjects regardless of similarity to a model bundle. HCP inter-subject variability was evaluated by generating voxel-wise heat-maps of the summed binary masks of each bundle.
Third, to test the variability of FWT in an independent dataset with different acquisition parameters and sampling schemes, we applied it to the MASSIVE dataset . Model-based wDSC bundle similarity scores, using the HCP-template bundles as references, were then calculated for every bundle from both the HCP and MASSIVE datasets. We computed wDSC descriptive statistics and generated violin plots for evaluation. This provided a framework for the comparison of the outputs from both datasets. Finally, we used single-rater agreement ICC scores to evaluate overall inter-session agreement.

Imaging data
2.2.1. HCP test-retest data-The test-retest HCP dataset (Van Essen et al., 2012) (https://db.humanconnectome.org/data/projects/HCP_Retest) consists of 2 scans acquired at 1 -11 months apart using the same scanning protocol on the same 3-Tesla Siemens Skyra scanner (Siemens Healthineers, Erlangen, Germany) using a 32-channel phased array receive head coil. The HCP T1-weighted images were acquired using a 3D Magnetization Prepared Rapid Acquisition Gradient Echo (MPRAGE) pulse sequence with 0.7 mm isotropic voxels. The diffusion data was acquired with 90 directions per shell with b-values (1000, 2000, and 3000 s/mm 2 ), and 1.25 mm isotropic voxels. We used the preprocessed (Andersson et al., 2003a;Sotiropoulos, 2016, 2015;Glasser et al., 2013;Jenkinson et al., 2002) imaging data of 20 random subjects (10 male / 10 female) resulting in 40 scans (see supplementary Table 1 for details). 25 years old) using various b-values and diffusion sampling schemes. Imaging data was acquired on a 3-Tesla Philips Achieva scanner (Philips Healthcare, Best, The Netherlands) with an 8-channel phased array receive head coil. MASSIVE 3D T1-weighted images were acquired with 1 mm isotropic voxels, while dMRI images were acquired with 2.5 mm isotropic voxels and multiple shells with (0-9000 s/mm 2 ) b-values. We used the (0-4000 s/ mm 2 ) b-shell data with the following b-values in s/mm 2 and number of diffusion-weighting gradient directions respectively (b500 -125, b1000 -250, b2000 -250, b3000 -250, b4000 -300). Five different sessions were generated using the preprocessed dMRI data acquired with an anterior-posterior (AP) phase-encoding axis and a negative gradient polarity to create 5 different sessions. Two sessions had multi-shell and three had single shell data with varying b-values, numbers of diffusion-weighted volumes and interleaved b0 volumes.
Additionally, we used corresponding reversed-phase encoded b0 images for Echo-Planar Imaging (EPI) distortion correction in FSL (Jenkinson et al., 2012). We used the MASSIVE data to investigate FWT reproducibility in the same subject with different b-values, number of diffusion directions and b-shells Table 1. lists the b-values and total number of volumes of all dMRI data used.

HCP test-retest group template creation-
We created group averaged T1weighted images and WM fODF maps from each individual HCP dataset in MRTrix3 (Tournier et al., 2019), referred to throughout this work as the HCP-template for simplicity. The resulting averaged T1-weighted image was parcellated using FreeSurfer (Fischl, 2012) and MSBP (Tourbier et al., 2019).

Virtual white matter dissection protocol-
We reconstructed the following 68 WM fiber bundles using FWT based on anatomical and neuroimaging literature definitions, 15 bilateral association bundles: The arcuate fasciculus, cingulum, fornix, frontal aslant tract, Inferior fronto-occipital fasciculus, inferior and middle longitudinal fasciculi, whole superior longitudinal fasciculus and its subcomponents, uncinate fasciculus, vertical occipital fasciculus. 8 commissural bundles: the anterior commissure, and corpus callosum in 7 segments. 12 bilateral projection bundles: the medial lemniscus, optic pathway, whole pyramidal tract and its subcomponents, and thalamic radiations. 3 bilateral cerebellar bundles: Dentato-rubro-thalamic tract, inferior and middle cerebellar peduncles. Definitions for each bundle are described in the results section and the inclusion/exclusion VOIs are provided in the supplementary methods. Fig. 1 provides a schematic description of FWT. A detailed technical description including tracking and optimization parameters is provided in supplementary methods.

Creation of inclusion and exclusion
VOIs-The first part of the FWT workflow automatically generates the inclusion and exclusion VOIs using the outputs of FS (Fischl, 2012) recon-all (FreeSurferWiki, 2020), MSBP (Tourbier et al., 2019), along with the UK BioBank (UKBB) volumetric atlas of fiber bundles , the spatially unbiased atlas template of the cerebellum and brainstem (Diedrichsen, 2006;Diedrichsen et al., 2011Diedrichsen et al., , 2009Diedrichsen and Zotow, 2015) (SUIT), the Neuroimaging and surgical technologies (NIST) Parkinson's disease histological atlas (Xiao et al., 2017(Xiao et al., , 2015(Xiao et al., , 2012) (PD25), and manually defined VOIs for the anterior and posterior commissures. The resulting anatomical VOIs are warped to diffusion space using ANTs (Avants et al., 2011), then combined to form bundle specific inclusion and exclusion VOIs. In our analysis, these VOIs were applied to each individual dataset (i.e., the two sessions of each subject in the HCP test-retest dataset, and the different scan sessions of the preprocessed MASSIVE single subject dataset). A similar workflow designed for group template T1 images was used on the HCP-template data.

Individual tractography and creation of FWT-HCP template bundles-
The second part of the FWT workflow provides the choice of two different whole brain WM virtual dissection approaches, both of which rely on the inclusion and exclusion VOIs created by the first part of FWT. The first approach employs bundle-specific seeding followed by streamline filtering and smoothing. The second generates a whole brain tractogram with 10 million streamlines by default followed by streamline dissection, then streamline filtering and smoothing. An equivalent workflow for fODF group averaged maps is also available generating whole brain tractograms with 20 million streamlines by default.
For this work we relied on bundle specific tractography to generate the 68 HCP-template bundles, and whole brain tractography followed by bundle dissection for all individual datasets. All template tractograms underwent visual quality assurance prior to further use. A tractogram is considered to have failed if less than 10 streamlines are generated initially, or if less than 10 streamlines are found after the first filtering step. We only included the successfully generated final FWT output tractograms in the current analysis.
2.3.7. Tractogram reproducibility measures-Tractograms can be used to generate segmentation maps that can be compared using well-known similarity/dissimilarity measures such as Hausdorff distance, overlap measures, DSC, etc. As streamline specific variants of these similarity measures have been shown to underestimate tractogram similarity  we used voxel-based weighted-dice similarity (wDSC) scores to simultaneously assess overlap and streamline density agreement per voxel in our inter-subject and intra-subject variability analyses (Cousineau et al., 2017). Additional similarity measures and their definitions, e.g., DSC, density correlation , volume overlap and overreach (Maier-Hein et al., 2017), and bundle adjacency (Garyfallidis et al., 2012) are provided in the supplementary material along with voxel-wise cumulative maps for each bundle to demonstrate inter-subject variability.

Qualitative bundle descriptions, template output and sample variability
Detailed qualitative descriptions of the virtual dissection protocol and demonstrative figures of bundles reconstructed from the HCP-template are provided below. Bundles are grouped by type (Meynert, 1888), i.e., association, commissural, and projection bundles, and cerebellar bundles are grouped separately. Radar-plots show the wDSC scores for each bundle, and wDSC summary statistics for each bundle group are included in supplementary table 3. Overall single-score agreement ICC was 0.713 (upper-bound = 0.789, lower-bound = 0.637, P < 0.05).

Association fiber bundles
3.1.1.1. Arcuate fasciculus (AF).: The AF is a major component of the dorsal language stream. It constitutes a perisylvian fronto-temporal pathway consisting of a long 'direct' association fiber system and two shorter 'indirect' bundles, see Fig. 2.
The direct bundle connects the ventral precentral, posterior inferior and middle frontal gyri with the middle and superior temporal gyri, and the two shorter bundles: (1) an anterior network connecting the supra-marginal and superior temporal gyri with the precentral gyrus, and (2) a posterior network connecting the posterior middle temporal gyrus with the angular gyrus (Bain et al., 2019;Bernard et al., 2019;Catani, 2006;Chen et al., 2015;Fernández-Miranda et al., 2015;Wakana et al., 2007;Wang et al., 2016;Yeh et al., 2018).
3.1.1.2. Cingulum (CG/Cing).: The cingulum bundle (CG) is the principal WM tract of the cingulate gyrus and the limbic system and is involved in a diverse range of functions spanning emotional, behavioral and sensorimotor control, mnemonic processing, nociception and executive function. Broadly speaking, the CG is a bidirectional fiber system that encircles the corpus callosum lateral to the cingulate gyrus, extending from the frontal lobes to the WM of the ventral temporal lobe with numerous lateral projections joining and leaving the bundle along its path (Jones et al., 2013;Metzler-Baddeley et al., 2012;Pascalau et al., 2018;Wu et al., 2016b), see Fig. 3.
The complexity of CG connections is inadequately characterized by dMRI, which tends to reconstruct either a single unilateral bundle extending from the subgenual cingulate to the temporal lobe, or two subdivisions encompassing the subgenual-dorsal, and dorsal-temporal segments. However, primate tracer studies suggest the CG could be subdivided into three or four regions based on differences in anatomical connectivity and neurotransmitter profile (Heilbronner and Haber, 2014). Here we reconstruct two sub-divisions in keeping with previous dMRI literature.

Cingulate Cingulum (CCing).:
This represents the dorsal (cingulate) component of the CG and is reconstructed by tracking the streamlines between the rostral and caudal anterior cingulate cortex and the posterior cingulate, isthmic posterior cingulate cortex and the precuneus. It is immediately superior to the body of the corpus callosum. It has a low curvature C shape, tipped on its open end with anterior/ventral extension to the subgenual cortex and ends posteriorly behind the splenium of the corpus callosum.

Temporal Cingulum (TCing).:
This represents the ventral component of the CG and is reconstructed by tracking the streamlines between the hippocampus and the isthmic posterior cingulate cortex and precuneus posteriorly. It is located in the medial temporal lobe, ascending from the anterior medial temporal lobe to the midline parietal region, behind the splenium of the corpus callosum.

Fornix (FX).:
The fornix is the major fiber pathway associated with the hippocampus, and comprises predominantly efferent fibers connecting the hippocampus with the prefrontal cortex, the anterior thalamic nuclei, the mammillary bodies, the ventral striatum, and the basal forebrain. Initially formed by the alveus and fimbria, WM of the bilateral hippocampus coalesce as the fornix crus and body before diverging again into the pre-commissural and post-commissural fornix columns which derive their name from their position relative to the anterior commissure. Recent research suggests fibers within the fornix are arranged topographically reflecting functional anterior-posterior gradient along the long axis of the hippocampus, with laterally located fibers arising from the anterior hippocampus and medially located fibers originating in the posterior hippocampus (Christiansen et al., 2017;Liacu et al., 2012;Pascalau et al., 2018;Strange et al., 2014), see Fig. 4. Here we reconstruct the fornix as a single lateralized bundle for each side. Reproducibility was poor and reconstruction failed in 19 datasets.

Frontal aslant tract (FAT):
The FAT is a recently described association bundle connecting the inferior and superior frontal lobe, and is involved in speech initiation, verbal fluency and executive function/inhibitory control. More specifically it connects the pars opercularis and pars triangularis in the inferior frontal gyrus (IFG) with the presupplementary motor area (pre-SMA), supplementary motor area (SMA), and the anterior cingulate cortex (Dick et al., 2019;La Corte et al., 2021;Pascual-Diaz et al., 2020), see Fig. 5.

Inferior fronto-occipital fasciculus (IFOF).:
The IFOF is a large, long-range association fiber bundle connecting the occipital and temporal lobes to the frontal lobes, specifically the lingual, posterior fusiform, cuneus and polar occipital cortex, with the inferior frontal gyrus, medial fronto-orbital region and frontal pole, see Fig. 6.
Notably, it narrows at the level of the extreme capsule. Though theoretically distinct from other temporal lobe association pathways, the IFOF runs in close proximity to the middle longitudinal fasciculus (MdLF), inferior longitudinal fasciculus (ILF) and uncinate fasciculus (UF), which may pose issues for some tracking algorithms (Forkel et al., 2014;Wu et al., 2016a). Considerable debate prevails regarding its exact functions however it is believed to serve the ventral visual and language streams along with the inferior longitudinal fasciculus (ILF) and uncinate fasciculus (UF) (Caverzasi et al., 2016). The right IFOF may be related to facial recognition functions and semantic visual stream processes (Herbet et al., 2018), while the left IFOF is related to semantic language functions (Almairac et al., 2015).

Inferior longitudinal fasciculus (ILF).:
The ILF is a large association tract connecting the occipital and temporal lobes and may play an important functional role in visual memory and emotional processing. It lies in direct contact with or close proximity to several bundles, including the UF, IFOF, AF, optic radiations (ORs) and tapetal fibers of the corpus callosum (CC). Whilst commonly depicted as a single bundle, the ILF may have up to four morphological subdivisions which reflect its occipital termination points, and include lingual, cuneate, fusiform and dorso-lateral occipital subcomponents. However, given the lack of consensus about the existence and functional significance of these subdivisions, here we construct the ILF as a single bundle (Herbet et al., 2018;Latini et al., 2017;, see Fig. 7.

Middle longitudinal fasciculus (MdLF).:
The MdLF is a large association bundle which is hypothesized to play a role in language, visual and auditory processing (De Witt Hamer et al., 2010;Dick and Tremblay, 2012;Kalyvas et al., 2020). Broadly speaking it is thought to connect the temporal pole, superior temporal gyrus, angular gyrus, superior parietal lobule and precuneus forming a distinct association bundle that runs medial to the AF and lateral to the IFOF (Makris et al., 2017(Makris et al., , 2009Maldonado et al., 2013;Menjot de Champfleur et al., 2013;Seltzer and Pandya, 1984;Wang et al., 2013), see Fig. 8.

Superior longitudinal fasciculus (SLF -I, II, and III).:
The SLF is a parietooccipital association fiber system, located in the dorso-lateral aspect of the cerebrum and generally believed to comprise three or four components Wang et al., 2016). Broadly speaking it connects the frontal with the occipital, parietal and temporal lobes, see Fig. 9. The SLF is related to language (Madhavan et al., 2014) visuo-spatial (Hong et al., 2019), and meta-cognitive (Zheng et al., 2020) functions. When tracked as a whole the SLF is highly reproducible, but may underestimate the frontal extent of the ventral division of SLF II.

SLF-I.:
The SLF-I connects the superior parietal lobule and precuneus with posterior superior frontal cortical areas. It is a rather short association bundle that is located above the level of the cingulate cortex (Wang et al., 2016). This was generated for all datasets on both sides but was the least reproducible of the SLF segments, most likely due to its small volume and low streamlines count. 1.1.8.2. SLF-II.: SLF-II is located more infero-lateral to SLF-I and originates in the anterior intraparietal sulcus and the angular gyrus, terminating in the posterior regions of the superior and middle frontal gyri. SLF-II can be subdivided into two distinct subcomponents (Barbeau et al., 2020), the first is the dorsal component (SLF-IId), which connects the supramarginal and inferior parietal cortices to the dorsal middle frontal gyrus. The second is the ventral component (SLF-IIv), which is longer, more ventral and connects the supramarginal and inferior parietal cortices to the rostral middle frontal gyrus.

SLF-III.:
SLF-III connects the intraparietal sulcus and inferior parietal lobule to the inferior frontal gyrus (Barbeau et al., 2020). It is immediately superior and medial to the arcuate fasciculus partially overlapping with the horizontal fibers of the Arcuate fasciculus.

Uncinate fasciculus (UF).:
The UF is typically characterized as a hook shaped bidirectional association fiber bundle linking the ventral, medial and orbital frontal lobes and rostral temporal lobes. Whilst DTI typically reconstructs a short hook-shaped fasciculus (Kurki et al., 2013), higher order dMRI reconstructions have revealed a more extensive fiber system which expands into a fan shaped trajectory in the frontal lobe. Such higher order models (e.g., CSD) in combination with microdissection suggest that the UF may be further subdivided into 5 subcomponents (Hau et al., 2017). However, for the present reconstruction, we opted to reconstruct the UF as one bundle for practicality, see Fig. 10.

Vertical occipital fasciculus (VOF).:
The VOF is described as a short slanted vertical association bundle in the lateral aspect of the occipital lobe. It connects the superior part of the occipital lobe and adjacent cortex of the occipito-parietal sulcus to the inferior aspect of the occipital lobe and adjacent occipito-temporal cortical areas (Jitsuishi et al., 2020;Schurr et al., 2019;Yeatman et al., 2014), see Fig. 11.

Anterior commissure (AC).:
The anterior commissure crosses the midline anterior to the pre-commissural columns of the fornix, above the basal forebrain and below the medial and ventral portion of the anterior limb of the internal capsule. It has two main parts, the first, anterior division which includes the olfactory decussation and the second, largest division which connects the temporal lobes, occipital lobes, and pre-dominantly the bilateral amygdalae (van Meer et al., 2016). During its course the AC intersects with the UF, ILF, sagittal stratum and optic radiations, which complicates tracking the true extent of the structure (Çavdar et al., 2021;Choi et al., 2011;Kikinis et al., 2015;Peltier et al., 2011;van Meer et al., 2016;Wilde et al., 2006), see Fig. 12.

Corpus callosum (CC).:
The corpus callosum is the largest fiber bundle in the human brain and connects the right and left cerebral hemispheres, hence its primary function is interhemispheric information transfer and integration. Traditionally the CC is divided into several subcomponents: the rostrum, genu, body, isthmus, splenium and tapetum. Here we reconstructed the CC in 7 segments comprising the prefrontal, premotor, motor, sensory, parietal, occipital and temporal callosal fibers, see Figs. 13 and 14.
It is common in DTI based literature to find three subdivisions: the genu, forming the forceps minor and connecting left and right prefrontal and anterior cingulate cortices; the callosal body, and the splenium, which forms the forceps major and connects left and right posterior parietal, medial occipital and medial temporal cortices. Improvements in fibertracking algorithms and integration with functional data have led to alternative parcellation strategies (Catani, 2006;Fabri et al., 2014;Nazem-Zadeh et al., 2012;Ohoshi et al., 2019Ohoshi et al., , 2019Phillips et al., 2013).

Prefrontal corpus callosum (CC PreF).:
This was reconstructed for all datasets using the bilateral prefrontal cortices and the genu and anterior third of the body of the corpus callosum.

Premotor and supplementary motor corpus callosum (CC PMC and SMA).:
This was reconstructed for all datasets using the caudal middle frontal gyri, the supplementary motor areas bilaterally and the central body of the corpus callosum.

Motor corpus callosum (CC motor).:
This segment connects the primary motor cortices (precentral gyri) of both hemispheres via the posterior third of the body of the corpus callosum. The motor CC was generated for all datasets.

Sensory corpus callosum (CC sensory).:
This was reconstructed using the primary sensory cortices (postcentral gyri) on both sides via the posterior third of the body and splenium of the corpus callosum in the midline. The sensory CC was reconstructed for all except 4 datasets.

Parietal corpus callosum (CC parietal).:
This was reconstructed for all data using the whole parietal lobes bilaterally and the splenium of the corpus callosum.

Occipital corpus callosum (CC occipital).:
This was reconstructed for all data using the occipital lobes on both sides and the splenium of the corpus callosum.

Temporal corpus callosum (CC temporal).:
This was reconstructed using the lateral aspect of the temporal lobes bilaterally and the splenium of the corpus callosum in the midline. The temporal CC reconstruction failed in 1 dataset.

Medial lemniscus (ML).:
The primary dorsal ascending tracts originate in the dorsal columns of the spinal cord, ascend through the medulla oblongata to the thalamic ventral-posterior medial (VPM) and lateral (VPL) nuclei and via the internal capsule to the primary somatosensory cortex in the postcentral gyrus (Jang and Seo, 2015;Peng et al., 2019a), see Fig. 15. The ML was reconstructed using the UKBB-derived posterior brainstem VOIs, the VPM and VPL thalamic nuclei, and primary sensory cortices.

Optic radiation (OR).:
The ORs connect the lateral geniculate nucleus of the thalamus to the primary visual cortex within the occipital lobe, see Fig. 16. They are critical for visual processing and present important challenges during surgical procedures involving the temporal lobe. Notably, Meyer's loop presents a particular challenge for tractography (Mehra and Moshirfar, 2021). This portion of the OR, which exhibits significant interindividual variability, is characterized by a sharp anterior projection in the anterior temporal lobe that then bends posteriorly to join the sagittal stratum. Moreover, within the anterior loop, OR fibers may intersect with fibers of the anterior commissure, ILF and tapetum, increasing the likelihood of spurious streamlines (Chamberland et al., 2018;Goga and Türe, 2015;Lim et al., 2015;Martínez-Heras et al., 2015). We generated two representations of the optic radiations, first using a classical anatomical definition with only the pericalcarine cortex as a cortical inclusion VOI, and second using the entire occipital lobe as a cortical inclusion VOI.

Optic tract (OT).:
The OT is formed of decussating axons from the contra-lateral optic nerve as well as non-decussating fibers from the ipsilateral optic nerve. It extends from the optic chiasm to the lateral geniculate nucleus of the thalamus (Hofer et al., 2010;Mehra and Moshirfar, 2021;Wu et al., 2012).

Pyramidal tract (PyT).:
The pyramidal tract fibers carry motor impulses from the cerebral cortex to the spinal cord through the brainstem. Here we reconstruct the Pyramidal tract (PyT) as a whole including the primary motor and sensory cortices along with the premotor, supplementary motor and parietal proprioceptive cortices superiorly and the whole brain stem inferiorly. We also provide dissections of this bundle based on specific sub-systems of the sensory-motor network (Chenot et al., 2019;Jeong et al., 2013;Peng et al., 2019b;Zhang et al., 2010), see Fig. 17.

Primary sensory-motor pyramidal tract (Corticospinal tract) (CST).:
The CST descends predominantly from the primary motor areas of the precentral gyrus and somatosensory areas of the postcentral gyrus through the corona radiata, posterior half of the posterior limb of the internal capsule, and cerebral peduncles to the rostral brainstem where it forms the medullary pyramids before crossing the pyramidal decussation on its further descent through the spinal cord. We provide two reconstructions of this bundle, one using the primary motor and sensory cortices (CST) and another using only the primary motor cortices, excluding the postcentral gyri (M1_CST). 1.3.4.2. Premotor pyramidal tract (PyT PMC).: This represents the pyramidal tract fibers originating in the dorsal premotor cortex (caudal middle frontal gyrus) and descending to the brainstem via the anterior half of the posterior limb of the internal capsule. This bundle had the lowest number of streamlines and was the least reproducible of the pyramidal tract components. 1.3.4.3. Supplementary motor pyramidal tract (PyT SMA).: This represents the pyramidal tract fibers originating in the supplementary motor area and descending to the brainstem via the centrum semiovale and the anterior half of the posterior limb of the internal capsule.

Thalamic radiations (TRs).:
The thalamic radiations are a group of 4-5 thalamocortical projection bundles comprising both afferent and efferent fibers. Collectively, the fibers of the thalamic radiations can be described as fanning out towards the cortex and banding up as they approach the thalami. As they run through the internal capsule and corona radiata they intersect with several other projection bundles, which in addition to the complex nature of thalamic microanatomy, can complicate reconstruction (Behrens et al., 2003;Kakou et  The posterior thalamic radiation is represented by the extended version of the optic radiation (OR OL) including the entire occipital cortex. All thalamic radiations were reconstructed using the whole thalamus as an inclusion VOI and excluding non-contributing thalamic nuclei, e.g., the VPM and VPL nuclei were excluded for the anterior thalamic radiation.

Anterior thalamic radiation (ATR).:
The ATR connects the dorso-medial (DM), dorso-lateral (DL) and anterior thalamic (ATN) nuclei to the prefrontal cortex and is thought to be involved in executive functions and complex planning (Niida et al., 2018). The ATR was reconstructed using the thalamus excluding all nuclei but the DM, DL, and ATN, and the ipsilateral medial and rostral prefrontal cortices.

Superior thalamic radiation (STR).:
The STR connects the ventral thalamic nuclei to the motor and sensory cortices via the superior thalamic peduncle, posterior limb of internal capsule, and the corona radiata, conducting cerebellar and basal ganglia input to the motor cortex (Bosch-Bouju et al., 2013). The STR was reconstructed using the thalamus, excluding the dorso-medial, dorso-lateral and pulvinar nuclei, and the ipsilateral primary motor, supplementary motor and dorsal premotor cortices (Younes et al., 2019).

Parietal thalamic radiation (PaTR).:
This was reconstructed using the entire parietal lobe as a cortical include and the thalamus excluding all except the VPM, and VPL nuclei. This resulted in a bundle of streamlines connecting the thalamic VPM and VPL predominantly to the primary sensory cortex.

Posterior thalamic radiation (PoTR/OR OL).:
This bundle was considered to be the same as the extended optic radiation using the entire occipital lobe, and was not reconstructed separately but will be added in future versions.

Dentato-rubro-thalamic tract (DRTT).:
The DRTT is the main efferent pathway from the cerebellum to the cerebral cortex and an important subdivision of the SCP in the context of deep brain stimulation-based neurosurgery. It originates in the dentate nucleus of the cerebellum, ascends through the brain stem where the majority of its constituent axons cross the midline at the superior cerebellar decussation (SCP) to synapse with the contra-lateral red nucleus (RN) in the midbrain (Coenen et al., 2020;Kwon et al., 2011;Mollink et al., 2016;Nowacki et al., 2018), see Fig. 19. From here the DRTT continues to the ventrolateral (VL) and ventromedial (VM) nuclei of the thalamus. Technically the DRTT terminates in the thalamus; however, for functional completeness we extend its trajectory to terminate in the primary motor cortex.
Whilst an ipsilateral component of the DRTT has been described in the literature, for the current purpose we have limited the reconstruction to the classical/contralateral DRTT. The left DRTT reconstruction failed in 4 datasets, while the right failed in 9 datasets. The overall DRTT wDSC showed a median of 0.498, minimum of 0.078 and Max.A.I.D. of 0.203. Additionally, we generated template DRTT bundles with 2,000 and 50,000 streamlines to illustrate the influence of the number streamline seeding attempts on tracking outcome and subsequent template-based reproducibility (supplementary Fig. 1).

Inferior cerebellar peduncle (ICP).:
The ICP consists mainly of afferent sensory fibers projecting from the spinal cord to the cerebellum. The ICPs are functionally involved in the maintenance of balance and posture via the integration of proprioceptive sensory and motor functions (Voogd, 2004). Microdissection studies reveal that it consists of four afferent bundles and 1 efferent bundle (Lingford-Hughes and Kalk, 2012), however here we reconstruct it as a single lateralized bundle using the dentate, fastigial and interposed nuclei of the cerebellum, and the ipsilateral medulla oblongata (Salamon et al., 2007;van Baarsen et al., 2016), see Fig. 19.

Middle cerebellar peduncle (MCP).:
The MCPs are large paired bundles that connect the brainstem to the cerebellum on both sides. They are often reconstructed in neuroimaging studies as a single commissural bundle connecting both cerebellar hemispheres (Leitner et al., 2015). The MCPs are thought to play a role in the modulation of skilled manual motor functions (Lingford-Hughes and Kalk, 2012) and have been shown to consist of three sub-fascicles (superior, inferior and deep) by microdissection studies (Jhaveri et al., 2018). We reconstructed the MCP as a whole for each side using the contra-lateral pons and cerebellar cortex. The resulting streamlines cross the midline in the pons, in line with its neuroanatomical definition. (Re et al., 2017;van Baarsen et al., 2016;Voogd, 2004), see Fig. 19.

HCP test-retest data results
In total we successfully reconstructed 2672 out of the 2720 attempted bundle reconstructions. The following bundles could not be reconstructed in all datasets: the DRTT (15 failures), premotor pyramidal tract (6 failures), optic tract (3 failures), fornix (19 failures), anterior commissure (2 failures), sensory CC (4 failures), and temporal CC (1 failure). Results of inter-session pair-wise similarity analysis per bundle for the different HCP test-retest subjects are shown as violin plots in Fig. 20.
Voxel-wise binary union maps per bundle show an expected pattern of maximum agreement in bundle cores and maximum variability around the periphery, supplementary Figs 2. and 3. Fig. 21 shows the wDSC scores derived from comparing the MASSIVE bundles to the HCPtemplate bundles superimposed on the results of the HCP test-retest dataset for reference. The single score agreement ICC was 0.917 (upper-bound = 0.944, lower-bound = 0.884, p < 0.05).

MASSIVE data results
In total we successfully reconstructed 316 out of 340 bundles (68 for each of the 5 MASSIVE datasets). 63 out of the 68 bundles were generated for every dataset, while the anterior commissure was only generated in 1 dataset, and the DRTT and fornix were not generated in any dataset.

Discussion
This work aims to facilitate the use of advanced tractography methods within the clinical research community by providing both an educational reference and the opensource tractography pipeline FWT for researchers and clinicians wishing to apply CSD tractography to typical HARDI datasets. We demonstrate the variability of the reconstructed bundles within and between subjects using two test-retest datasets and offer solutions to improve tracking results for the least reproducible bundles.

Variability of virtual dissections
Our results demonstrate that inter-bundle variability was much greater than inter-subject variability for the same bundles i.e., there was a wide range of similarity in the extent to which different bundles were reconstructed irrespective of the dataset used. The AF and PyT for example were reconstructed in largely the same way in all datasets, whereas the fornix and DRTT showed high variability in density and spatial extent and were sometimes not reconstructed at all. In line with this, inter-subject agreement was highest for bundles with high streamline densities e.g., PyT, AF, CC and lowest for bundles with low streamline densities e.g., fornix, and DRTT.
Inter-session and intra-subject agreement were high for all bundles i.e., bundles were reconstructed largely to the same extent, (or failed to be reconstructed) in both sessions of the HCP test-retest dataset. In the test-retest data from a single-subject with varying b-values and shell schemes (MASSIVE) changing the b-value(s), number of gradient directions or shell scheme did not consistently alter bundle similarity to the model.
In agreement with other studies our results show considerable variability, influenced not only by data acquisition parameters such as scanner manufacturer, b-value, diffusion sampling scheme, reconstruction model (Schilling et al., 2021a), but also by the virtual dissection protocol itself (Schilling et al., 2021b). In the context of our study, despite using exactly the same WM dissection protocol for the HCP-template as for the individual datasets, none of the tracked bundles failed in the template data, while some failed in the individual datasets. This highlights the importance of data quality, as group-averaged data tends to have a higher signal-to-noise ratio (SNR) and increased blurriness in fine cortical structures due to inter-subject variability. These results demonstrate that even in a unified fully automated WM dissection protocol considerable inter-bundle and inter-subject variability remain. In contrast to previous work, we found that data acquisition parameters and sampling schemes had only a limited impact on bundle variability. This is likely due to the high number of volumes used per session, which would result in a higher SNR compared to the typical single-shell acquisition datasets used in clinical settings.
In agreement with other studies (Bonilha et al., 2015;Cousineau et al., 2017;Gu et al., 2017;Zhang et al., 2019), we found that larger bundles, particularly if they have both a larger volume and a higher streamline count, tend to be more similar/reproducible regardless of differences in scanning parameters or even normal inter-subject anatomical variations. It follows that automated procedures are most appropriate for fasciculi with these characteristics. For smaller, less dense or more narrow bundles such as the fornix, DRTT (and relatedly, the superior cerebellar peduncle) or specific subdivisions of larger bundles (e.g., SLF-I, temporal corpus callosum) we recommend using either data with sub-millimeter spatial resolution, bundle-specific manual dissection or using the provided template bundles in an automated RecoBundles (Garyfallidis et al., 2018) based workflow rather than the VOI-based approach. Of note, we also provide symmetrical versions of the HCP-template bundles along with aligned reference anatomical images. In line with the use of template bundles, FWT is compatible with RecoBundles (Garyfallidis et al., 2018) based pipelines that can overcome the need for structural parcellation, which can be particularly useful in clinical data to avoid false tractography results due to structural parcellation errors.

The benefits of using a CSD atlas in clinical studies
There are currently a number of diffusion-based tractography atlases and protocols available (Catani and Thiebaut de Schotten, 2008;Thiebaut de Schotten et al., 2011;Verhoeven et al., 2010;Wakana et al., 2004;Warrington et al., 2020;Wasserthal et al., 2018;Yendiki et al., 2011) some of which are based on HCP data Yeh et al., 2018). Here we aimed to primarily address the limitations of DTI based virtual dissection protocols, which predominate in clinical studies, either as part of comparative studies with healthy controls, or in a healthcare setting using data acquired in shorter scan times than are achievable in a research environment. We did not include formal comparisons to other HARDI/higher-order modelling-based atlases which precludes us from drawing conclusions about the relative benefits of FWT over other methods for virtual dissection in this context. It is our aim to complement other atlases based on detailed anatomical research and specialized acquisition schemes with a comprehensive description of CSD-based human WM anatomy which improves upon classical DTI representations in 3T data. Whilst the results presented here are based on high quality data in healthy subjects, we also tested our protocol on lower quality clinical data during the development phase of FWT, to ensure its translatability to clinical populations (see supplementary Fig. 9). A formal analysis relating to the application of FWT in patients with lesions is beyond the scope of this paper focusing on typical CSD anatomy and will be presented in a follow-up publication.

Technical considerations and limitations
This work relied on imaging data from only 20 individuals, however by using re-test scans yielding 40 datasets and the additional analysis using the MASSIVE dataset, our sample is sufficient to study test-retest reproducibility, as well as inter-subject variability. While our dissection protocol was literature-based and guided by several contemporary publications we did not include a direct comparison to any other atlas, nor did we apply formal specific criteria for assessing the validity of the definitions. Given the lack of consensus on anatomical definitions relating to WM fasciculi, it is possible that some readers will disagree with our interpretation, and that the FWT VOIs and reconstructions may need to be revised as the field evolves.
We attempted to keep rigid heuristic decisions (e.g., statistical thresholds) to a minimum especially those pertaining to streamlines filtering, where we utilized largely data-driven streamline filtering tools from Dipy (Garyfallidis et al., 2014) rather than relying on a more stringent selection of anatomical exclude VOIs. However, in order to cater to varying data quality, our choice of inclusion VOIs in some situations are influenced by considerations made for lower quality datasets than the ones utilized in this work. For example, for the AF we excluded any streamlines involving the precentral gyrus to avoid false streamlines tunneling through volume average voxels, which we have observed in lower quality singleshell data. Similarly, the AC is generated without the olfactory component, and we do not include subdivisions to the AF, CST, full posterior thalamic radiation, or superior cerebellar peduncles proper. We also do not include more esoteric bundles such as the superior-anterior fasciculus (David et al., 2019), temporo-insular fasciculus Radwan et al., 2019), and inferior (auditory) thalamic radiation (Maffei et al., 2019a(Maffei et al., , 2019b, or those with a high chance of failure due to their geometry, such as the posterior commissure, or cranial nerves. FWT tractography utilizes spherical-deconvolution informed filtering of tractograms (SIFT) (Smith et al., 2015(Smith et al., , 2013) based streamlines seeding and filtering for whole brain tractography, while further improvement can be expected if an optimized whole brain tractography method is utilized, e.g., global tractography (Christiaens et al., 2015) or particle-filtering tractography (Girard et al., 2014). FWT does not use a specific required number of output streamlines during bundle segmentation from whole-brain tractograms. However, this is necessary if the bundle specific approach is utilized. In this case the resulting bundles will be directly influenced by the required number of streamlines as this will influence the number of streamlines seeding attempts. This point is illustrated in supplementary Fig. 1, which shows that when the DRTT template bundles are generated with 2,000 and 50,000 required streamlines, the outcome from the 50,000 streamlines appears more complete as expected. Furthermore, when applied to the individual datasets, the bundle-specific approach with 20,000 required streamlines reconstructed more bundles successfully than the whole-brain tractography and bundle segmentation approach (73 vs. 67), and resulted in higher wDSC scores particularly when compared to the 50,000 streamlines DRTT template bundles. Tracking such challenging bundles may also benefit from anatomically-constrained tractography , which will be adopted in future versions of FWT.
The FWT pipeline is time-consuming, as FS recon-all alone requires at least 4-6 h and dMRI preprocessing can vary widely depending on acquisition parameters, and reconstruction method. E.g., whole brain tractography can take more than 4 h for a multishell multi-tissue CSD model with probabilistic tractography using second-order integration over distributions of fODFs (iFOD2) (Tournier et al., 2010) with fODF-driven dynamic seeding (Smith et al., 2015). Whilst this is offset by the flexibility of FWT allowing a more targeted approach using bundle-specific seeding and tractography in case of presurgical mapping, further development will focus on improving processing time. Additionally, future work will provide a clustering-based workflow independent of prior structural parcellation, the addition of more fiber bundles, as well as generating deterministic versions of all bundles. Finally, any justified and necessary changes to the inclusion/exclusion VOIs for any bundle can be easily implemented by changing the FWT workflows, which are openly provided, and a new atlas of bundles can be easily created by rerunning the template workflows on the HCP-group template or similar data.

Conclusion
The FWT pipeline can reconstruct 68 WM fasciculi and shows high inter-and intra-subject reproducibility. Dense bundles such as the PyT, AF and MCP yield the most reliable and least variable reconstructions. Higher resolution data or tract specific modification may be required for thinner bundles such as the fornix, DRTT, SLF I and premotor PyT. The FWT CSD atlas may be a useful reference and virtual dissection tool for applied neuroimaging students and clinical professionals wishing to use and understand the capabilities of CSD tractography.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material. Graphical representation of the 2 parts of FWT for the corticospinal tracts (CSTs), arcuate fasciculi (AF), and motor corpus callosum (CC). FWT part 1 (top) requires as input preprocessed diffusion data, and FS recon-all output, and MSBP output. The script generates all VOIs used for virtual dissection by combining various anatomical VOIs from different parcellation maps and atlases. FWT part 2 (bottom) generates all tractograms from preprocessed diffusion data and VOIs created by FWT part 1, this script provides two approaches to virtual dissection; (1) Bundle specific and (2) Whole brain tractography followed by segmentation. FS = FreeSurfer, MSBP = MultiScale Brain Parcellator, PD25 = NIST Parkinson's histological, JHU = John's Hopkins university, Juelich = Juelich university histological atlas, AC/PC = anterior commissure/posterior commissure, manually defined VOIs in template space, UKBB = UK Biobank, SUIT = spatially unbiased cerebellar atlas template, dMRI = diffusion magnetic resonance imaging, CSD = constrained spherical deconvolution, fODF = fiber orientation distribution function. $ = preprocessed, should include correction for motion, Eddy currents, EPI distortion, imaging noise and bias. * = Other models e.g., DTI with FACT.                  Table 1 Diffusion weightings and total number of volumes.