Quantitative susceptibility atlas construction in Montreal Neurological Institute space: towards histological-consistent iron-rich deep brain nucleus subregion identification

Iron-rich deep brain nuclei (DBN) of the human brain are involved in various motoric, emotional and cognitive brain functions. The abnormal iron alterations in the DBN are closely associated with multiple neurological and psychiatric diseases. Quantitative susceptibility mapping (QSM) provides the spatial distribution of the magnetic susceptibility of human brain tissues. Compared to traditional structural imaging, QSM provides superiority for imaging the iron-rich DBN owing to the susceptibility difference existing between brain tissues. In this study, we constructed a Montreal Neurological Institute (MNI) space unbiased QSM human brain atlas via group-wise registration from 100 healthy subjects aged 19–29 years. The atlas construction process was guided by hybrid images that were fused from multi-modal magnetic resonance images (MRI). We named it as Multi-modal-fused magnetic Susceptibility (MuSus-100) atlas. The high-quality susceptibility atlas provides extraordinary image contrast between iron-rich DBN with their surroundings. Parcellation maps of DBN and their subregions that are highly related to neurological and psychiatric pathology were then manually labeled based on the atlas set with the assistance of an image border-enhancement process. Especially, the bilateral thalamus was delineated into 64 detailed subregions referring to the Schaltenbrand–Wahren stereotactic atlas. To our best knowledge, the histological-consistent thalamic nucleus parcellation map is well defined for the first time in the MNI space. Compared with existing atlases that emphasizing DBN parcellation, the newly proposed atlas outperforms on the task of atlas-guided individual brain image DBN segmentation both in accuracy and robustness. Moreover, we applied the proposed DBN parcellation map to conduct detailed identification of the pathology-related iron content alterations in subcortical nuclei for Parkinson’s Disease (PD) patients. We envision that the MuSus-100 atlas can play a crucial role in improving the accuracy of DBN segmentation for the research of neurological and psychiatric disease progress and also be helpful for target planning in deep brain stimulation surgery.


Introduction
Three-dimensional standard space brain atlases are valuable tools as an anatomical reference for human brains (Horn and Kühn 2015;Fonov et al. 2011;Mazziotta et al. 2001). Atlases are widely used to measure neurological disorder related structural changes in patients to compare with the healthy population (Mazziotta et al. 1995;Xiao et al. 2012;Fonov et al. 2009;Butson et al. 2007). Brain atlases that emphasize subcortical Deep Brain Nuclei (DBN) structures can also be used to refine targets of neurological disease treatment, such as deep brain stimulation (DBS) (Merkl et al. 2016;Neumann et al. 2015;Eisenstein et al. 2014). DBN are closely related to the pathogenesis of neurodegenerative diseases, such as Parkinson's disease (PD), dystonia, or obsessive-compulsive disorder. Faithful identification and depiction of DBN based on neuroimaging techniques are crucial for in vivo brain structural analysis (Lenglet et al. 2012;Starr et al. 2002;Voges et al. 2002). However, at 3T, conventional T1-weighted (T1w) and T2-weighted (T2w) MR images have difficulties in clearly delineating the boundaries of DBN from surrounding tissues. As well, the image Chenyu He and Xiaojun Guan contributed equally to this work.
* Yuyao Zhang zhangyy8@shanghaitech.edu.cn Extended author information available on the last page of the article contrast for substructures within DBN is also homogeneous. This is mainly because tissue contrast between structures in T1w/T2w originates from the differences in signal relaxation rates of the water protons, which are rather weak in DBN. Thus, the existing T1w atlases may not be ideal for the depiction of DBN in the human brain.
There have been several atlases that were constructed with specific emphasis on the localization of brain subcortical nuclei. For example, the DBS intrinsic template atlas (DISTAL atlas) (Ewert et al. 2018) provides subcortical anatomy on ICBM 152 Nonlinear 2009b templates. In this work, three of the most common DBS targets (subthalamic nucleus (STN), internal part of the pallidum (GPi), and ventral intermediate nucleus of the thalamus (VIM)) are localized by co-registering a detailed histological atlas onto the standard T2w template in MNI space. However, the spatial normalization between the individual histology map and the standard space unbiased T2w template may induce registration deviation. Based on the construction of a high-quality T2w unbiased template, the California Institute of Technology (CIT168) atlas provides probabilistic anatomical labels for a bunch of subcortical nuclei that are related to the reward learning and decision making circuits in the human brain (Pauli et al. 2018). In the CerebrA atlas (Manera et al. 2019), the authors combined the ICBM 152 Nonlinear 2009b symmetric model (at a resolution of 1 ×1× 1 mm 3 ) with the Mindboggle-101 (Klein and Tourville 2012), for propagating the labels from Mindboggle-101 to the MNI space through T1w template with an additional manually correction on the subcortical labels. Lau et al. (2020) demonstrated that high-resolution (sub-millimetric) T1 relaxometry measurements at ultra-high field strength (7T) could be used to delineate zona incerta (ZI), a small nucleus inbetween STN and red nuclei (RN), from surrounding white matter structures. Using this approach, the authors successfully derived in vivo estimates of the size, shape, location, and tissue characteristics of ZI, confirming observations that were previously only estimated through histological evaluation. Notably, all the studies mentioned above benefited from the susceptibility-sensitive contrast (T2w or longitudinal T1 relaxometry measurements) that provided better image contrast between DBN and their surrounding tissues.
Recently, a novel post-processing method, referred to as quantitative susceptibility mapping (QSM) (Liu 2010;Schweser et al. 2011;Shmueli et al. 2009;Wharton et al. 2010) has been introduced. QSM provides the quantitative measures for iron content of human tissue composition (Bilgic et al. 2012;Langkammer et al. 2010). QSM is sensitive to the spatial variations of molecular or cellular components that exhibit different magnetic susceptibility properties. Due to its high sensitivity to iron and myelin, QSM provides clear tissue boundaries between iron-rich DBN and myelinated white matter. Thus, QSM enabled the identification of several substructures of the human brain that are in part indiscernible in other MRI contrasts (i.e., T2w, R2*). As indicated by several recent studies (Shmueli et al. 2009;Haacke et al. 2010;Schweser et al. 2011;Hanspach et al. 2017), common DBS targets, including substantia nigra pars reticulata (SNr), VIM, STN and the substructures of GP, are clearly distinguishable from surrounding tissues in QSM. Thalamus also attracted plenty of attention in QSM-related studies, such as in Schweser et al. (2018) and Hanspach et al. (2017). Both authors proposed to depict thalamic subnuclei directly based on the contrast of the group-wisely averaged QSM template. Deistung et al. (2013) indicated that using a high-resolution QSM image acquired at 7T, QSM exhibited superb contrast that allowed histological-consistent identification for substructures of the thalamus, midbrain, and basal ganglia.
There are a few susceptibility atlases proposed for facilitating DBN parcellation (Lim et al. 2013;Zhang et al. 2018;Li et al. 2019b). Lim et al. (2013) added a QSM template to the single-subject Eve atlas space, to extend the Eve atlas with an accurate "deep gray matter parcellation map". Limited by the low signal-to-noise ratio (SNR) of the single QSM image and the lack of consideration of heterogeneity among subjects, the parceled DBN could be biased and less representative. Zhang et al. (2018) proposed a longitudinal age-dependent QSM atlas, providing an efficient tool for segmenting brain subcortical structures from individual QSM images for a variety of age ranges. However, restricted by the image resolution (1 × 1 × 1 mm 3 ) of the 3T QSM template, the thalamus was parceled into only five subregions in each hemisphere. Li et al. (2019b) proposed a multi-atlas strategy for achieving a more reliable automated tissue susceptibility quantification framework in a series of typical DBN. However, to the best of our knowledge, an MNI space unbiased version, and a histological-consistent DBN subregion parcellation map are not available in these QSM templates.
For the standard MNI space, the most widely used atlases are Colin-27 atlas (Holmes et al. 1998) and ICBM 152 atlas series (Mazziotta et al. 2001;Fonov et al. 2009Fonov et al. , 2011. Colin-27 includes T1w/T2w atlases with a normal-resolution version (1 × 1 × 1 mm 3 ) and a high-resolution version (0.5 × 0.5 × 0.5 mm 3 ). In the ICBM 152 nonlinear 2009b version, T1w/T2w/PD and T2-relaxometry atlases are available with both symmetric and asymmetric versions. Especially, the two groups of atlases are both constructed from normative young adult populations: Colin-27 (a 28-years individual) and ICBM 152 (18.5-to-43.5-years individuals). This is partially because the MR image quality is sensitive to motion artifacts during image scanning. To minimize the effect of motion on image quality and thus further facilitate the atlas construction process, normative young adult subjects are adopted in adult brain atlas construction research works.
Towards this end, we proposed to construct an unbiased QSM atlas in the MNI space, using multi-modal MR images (including modalities of T1w, Gradient Echo (GRE) magnitude, and QSM) collected from one hundred young healthy individuals (age 19-29 years). In our strategy, T1w images provide clear GM-WM contrast and act as the guidance towards the MNI space; GRE magnitude is used for normalizing T1w and QSM spatial information; while QSM performs superb contrast for inside and outside DBN and highlights positions for iron-rich nuclei, and especially guides detailed DBN subregion identifications in the MNI space. The averaging process for spatially normalized individual MR images implemented in the atlas construction substantially enhanced the image SNR in the unbiased QSM template. We further applied the border-enhancement technique to emphasize the boundaries of DBN and then proposed a series of DBN delineations, including 32 thalamic subnuclei on each hemisphere of the thalamus (Thal), and another ten pairs of bilateral basal ganglia nuclei: putamen (Pu), nucleus accumbens (NAcc), caudate nucleus (Cau), ventral pallidum (VP), internal and external globus pallidus (GPi and GPe), red nucleus (RN), pars reticulata and pars compacta of substantia nigra (SNr and SNc) and subthalamic nucleus (STN). Especially, the depiction of thalamic subnuclei was conducted on a handcrafted parcellation map on an individual 7T QSM image following the Schaltenbrand-Wahren's stereotactic atlas (Schaltenbrand 1977), then carefully warped into the 3T individual space, and weightily fused into the MNI space. To our best knowledge, the detailed histologicalconsistent thalamic nucleus parcellation map is defined in the MNI space for the first time. To assess the feasibility of the MuSus-100 atlas, we additionally segmented six individual subjects using the proposed atlas and then calculated quantitative parameters by comparing segmentation results with ground truth (GT) manual annotations. The result shows that the MuSus-100 atlas outperforms existing atlases on atlas-guided individual brain image DBN segmentation, both in accuracy and robustness. We also applied the proposed atlas on the result of voxel-based analysis (VBA) PD QSM image analysis, to identify the pathology-related iron content difference in small-size DBN subregions between PD patients and healthy control subjects. Abnormalities in iron content are found in the bilateral Cau, Pu, GP, Thal, RN, left STN, and SN (both SNr and SNc). In the thalamus, the ventral and the central subregions show the most essential PD-related brain tissue susceptibility variation.

Methods
This section describes the detailed processes of MNI space unbiased QSM atlas construction and the parcellation of DBN, mainly including steps of participant recruitment, image acquisition, atlas construction, DBN parcellation, evaluation criterion, and applications of the proposed atlas.

Data for atlas construction
One hundred healthy volunteers (41 males, 59 females) were recruited for MR data collection. The volunteers were aged between 19 and 29 years, with the mean value of 24.03. The data were acquired on a 3T GE750 MRI scanner. Each individual underwent the scanning sessions of a T1-weighted (T1w) sequence and a multi-echo, spoiled-gradient recalled echo (SPGR) sequence. The scanning parameters of the T1w images: repetition time (TR) = 7.336 ms; echo time (TE) = 3.036 ms; flip angle = 11 • ; field of view (FOV) = 260×260 mm 2 ; matrix size = 256×256; slice thickness = 1.2 mm; 196 continuous sagittal slices. A three-dimensional SPGR sequence was utilized to obtain QSM images with the following parameters: TR = 33.7 ms, TE1/spacing/TE8 = 4.56/3.65/30.11 ms, flip angle = 20 • , FOV = 240×240 mm 2 , matrix size = 416×384, slice thickness = 2 mm, resolution = 0.47×0.47× 2 mm 3 . All the images were resampled to the same resolution of 1 ×1× 1 mm 3 through the operations in k-space. The raw phase was unwrapped using the Laplacianbased phase unwrapping, and the normalized background phase was removed by V-SHARP (Wu et al. 2012). The susceptibility maps were determined by the STAR-QSM algorithm (Wei et al. 2015). All scans were inspected to screen out the images with incomplete brain structure or with obvious motion artifacts.
All the procedures of the present study are in accordance with the Declaration of Helsinki and are approved by the Ethics Committee of the Second Affiliated Hospital, Zhejiang University School of Medicine. The participants were all informed of the study procedures and gave their written consent.

Atlas construction
The atlas construction mainly relies on registrations between individual images. To make the cortical and iron-rich subcortical regions both accurately aligned in registrations, we fused the individual T1w and QSM scans to create the hybrid images. And the atlas construction process was performed on the hybrid instead of the QSM individual scans, to achieve accurate structural alignment with sharp image contrast in both cortical and subcortical regions. Thus, the atlas construction process mainly includes four steps: individual T1w-QSM alignment, hybrid image generation, native space atlas construction, and MNI space projection. The pipeline is illustrated in Fig. 1, where the individual T1w-QSM alignment and hybrid image generation are shown in the pink fusion module (described in "Hybrid individual image construction"), and the native space atlas construction (described in "Native space atlas construction") along with the MNI space projection (described in "Towards MNI space projection") are noted as (a) and (b) parts in the figure.

Pairwise image registration
The pairwise image registration process takes an image pair as input, named moving image ( I M ) and fixed image ( I F ). The process aligns the moving image to the fixed image by performing proper transformation to the moving image. The alignment includes the coordinate space and the visual representation between images. The transformation is calculated by minimizing the dissimilarity between the transformed moving image and the fixed image. Thus, the registration process includes two steps, i.e., the dissimilarity optimization step and the transformation application step. The dissimilarity optimization step can be expressed as where L sim and L smooth , respectively, denote the dissimilarity metric and the regularization term for keeping the transformation smooth. With the yielding transformation M� →F , the image transformation step can be expressed as Combining the dissimilarity optimization and image transformation steps defined in Eqs. 1 and 2, the general registration process can be formulated as where the output of a pairwise registration includes the aligned moving image I M� →F and the transformation M� →F .
The pairwise registration process can be achieved by either linear or nonlinear transformations, formulated as linear registrations or nonlinear registrations. Usually, linear transformations are described via transformation matrices, while nonlinear ones are described via the deformation fields. To improve accuracy, nonlinear registrations typically require linear registration for initial alignment. We present the combination of linear and nonlinear transformations into a single function f reg (I F , I M ) as shown in Eq. 3.
In practice, we used affine linear registration as the initialization, followed by a symmetric diffeomorphic (SyN) (Avants et al. 2008) for the nonlinear registration process. Before the registration step, T1w images underwent bias field correction and skull removal processing. Bias field correction was conducted using N4 intensity normalization (Tustison et al. 2010), and skull removal was performed using Brain Extraction Tool (BET) (Smith 2002).

Hybrid individual image construction
In this section, we describe the motivation and the detailed process of constructing individual hybrid images. The process is marked as the pink fusion model in Fig. 1, where the detailed process is illustrated in the pink box located in the top right of the figure.

Motivations of constructing hybrid images
The primary process of the atlas construction is the group-wise registration. The accuracy of pairwise registrations affects the quality of the constructed atlas. High sharpness of the boundary of the cortical and subcortical areas would facilitate an accurate alignment of registration. We did not use individual T1w images for the group-wise registration because of its poor contrast in subcortical regions, leading to potential misalignment of DBN. Instead, we first fused the individual T1w and QSM images to obtain individual hybrid images. The hybrid images display sharp contrasts in both cortical and subcortical regions, suitable for the DBN-focused atlas construction.
Native space T1w and QSM image space alignment Each individual has multi-modal MR images of T1w, QSM, and GRE magnitude (Mag). Among these, QSM and GRE magnitude images are reconstructed from multi-echo GRE sequences, whereas T1w structural images are reconstructed from fast spin-echo sequences. Thus the T1w and the QSM images are not in the same coordinate space, while QSM and GRE magnitude images share the same coordinate space. GRE magnitude images have relatively sharp structural contrast, which is suitable for registering to the T1w images. Thus, we used GRE magnitude images as the connection between T1w and QSM images.
The process of individual T1w-QSM alignment is denoted as follows. For each individual i ( i ∈ [1, N] , N denotes the total number of the subjects), the GRE magnitude image I Mag i is aligned to the T1w image I T1w i using pairwise image registration, which is The pipeline of MNI space unbiased QSM atlas construction.
(a) Native space atlas construction. Individual hybrid images are fused from individual T1w and QSM images. The fusing module is shown in the right-top box (pink), indicating the pairwise image registration between individual T1w and QSM images for generating the hybrid image. Group-wise registration is then applied to individual hybrid images to generate a hybrid atlas and warp fields for each subject. The individual T1w and QSM images are projected to the corresponding native atlas space by applying the warp fields from the group-wise registration. The native space T1w and QSM atlases are constructed by averaging the individual T1 (green) and QSM (orange) images registered to the group mean in the corresponding native space, respectively. (b) Atlas projection into the MNI space. The constructed native space T1w atlas (green) is projected to the MNI space (dark green) via registering to the ICBM T1w atlas (black), yielding the transformation from native space to the MNI space. The native space QSM atlas (orange) is projected to the MNI space (brown) via applying the transformation 1 3 Then, the transformation Mag� →T1w i is applied to the individual QSM image I QSM i to transform it to the coordinate space of T1w image, expressed as Generally, the operation of multi-modal image alignment process, combining the registration process and the transformation process denoted in Eqs. 4 and 5, can be expressed as After the pairwise registration, hybrid images are generated by fusing T1w and QSM images. For each individual i, the fusion process is a weighted linear combination of the input image pairs (I T1w where is a scalar controlling the fusion extend. Increasing makes the result carries more visual feature from QSM, thus obtaining a stronger visual contrast of DBN, and vice versa. After multiple trials, = 12500 is empirically selected by visual inspection. At this setting, hybrid images well balance the contrast in cortical and subcortical regions.

Native space atlas construction
Native space atlas construction is a process of deriving an averaged template according to a dataset of MRI scans, shown in Fig. 1a. The averaged template is achieved by a group-wise registration, followed by a template update process to sharpen the constructed atlas. We performed the atlas construction process on hybrid images yielding the native hybrid atlas with individual transformations, and then we constructed the native T1w and QSM atlases by averaging the warped multi-modal individual images.
Unbiased group-wise registration is the main step of the atlas construction. The process takes a set of images as input, deriving an unbiased average template after several iterations. Initially, we chose the individual image that is the most similar to all the other images in the input dataset as the initial template T (t=1) . For each iteration, the algorithm aligns each individual image to the current averaged template T (t=m) (where (t = m) denotes the m-th iteration) constructed in the last iteration, yielding the transformation for each individual to the current template. The individual-template alignment is achieved by pairwise registrations denoted in "Pairwise image registration". Then, each original individual image I (t=0) i is updated to the current template T (t=m) at the m-th iteration by registering to the template as Finally, the template is updated by voxel-wise averaging of transformed individual images (t=m) as Generally, we mark the group-wise registration process at the m-th iteration as Post-processing template update Additional post-processing steps are applied to the average atlas to achieve a refined template update. We mark the process as , where (t=m) −1 represents the inverse transformation set at the m-th iteration. The detailed process includes: 1. Voxel-wise averaging of the nonlinear inverse warp fields, which can be considered as the level of brain structure variation among the population used for atlas construction. 2. We further multiplied a factor for controlling the image deformation level on the averaged inverse warp fields. The factor is a scalar defined in relation to the gradient step, which is a predefined hyperparameter controlling the image deformation level in SyN (Avants et al. 2008). The range of the scalar is [0, 0.25] and the default value of 0.2 was used. 3. Applying transformations to the averaged atlas, including applying the inverse linear transformations once and applying the averaged nonlinear inverse warp field four times.
Native space hybrid atlas construction By implementing the above two steps on the hybrid image set, we acquired the hybrid atlas, along with the transformations for each individual image to the constructed hybrid atlas space. In each iteration, the full process of native space hybrid atlas construction, including hybrid image generation, group-wise registration, and average template update, is summarized in Algorithm 1. The algorithm (excluding the hybrid image construction) was performed using antsMultivariateTem-plateConstruction2 utility implemented in the Advanced Normalization Tools (ANTs) registration toolkit (Avants et al. 2009).
Algorithm 1 Native space hybrid atlas construction Choose the individual image with the minimal intracluster distance 4: while m ∈ {1, · · · , M} do M denotes the total iteration step 5: Original scan warped to current template 7:T (t=m+1) ← avg(I hybrid; (t=m) ) 8: end while 10: return T, Φ 13: end procedure Native space T1w and QSM atlas construction The T1w and QSM atlases were created by applying the transformations from the atlas construction process to the corresponding individual images and averaging the warped individual scans, which is denoted as Note that only the averaging process is performed while the update process in atlas construction is not performed, since updating is a plug-in process during the atlas creation process.

Towards MNI space projection
For projecting our unbiased QSM atlas into the standard MNI space, we first performed pairwise registration between our T1w atlas and the modal-identical standard MNI space T1w atlas (ICBM Nonlinear Asymmetric 2009c, resolution of 1 ×1× 1 mm 3 , short for ICBM atlas) (Ewert et al. 2018). For ICBM T1w atlas, the skull was removed using the official brain tissue mask. Then, we warped our T1w atlas towards ICBM T1w atlas via pairwise registration. In this process, the deformation field was saved and then applied to the native space unbiased QSM atlas. The process is shown in Fig. 1b.

DBN border-enhancement and manual parcellation
For further facilitating DBN depiction using the proposed atlas, we applied border enhancement by applying a 3D volumetric Laplacian kernel to the constructed QSM atlas. The implemented Laplacian kernel is defined as The above kernel was applied to both QSM and T1w templates to yield QSM and T1w Laplacian maps. The image intensities of the Laplacian maps were scaled to [0,255]. For annotating DBN, the generated templates and Laplacian maps were both taken as references, whereas the QSM template and its Laplacian map were mainly referenced due to the clear contrast in iron-rich DBN in those references.
Two experienced radiologists participated in the delineation work, to eliminate subjectivity from single-rater delineation. Considering that a common initialization often leads to a more consistent manual delineation between different raters, a full-brain parcellation map AAL3 (Rolls et al. 2020) that contains DBN region annotations is first overlapped on the MNI space unbiased QSM template. AAL3 contains fewer DBN regions than we propose, the rest of the ROIs are defined based on QSM template contrast. The delineation work was carried out by the two raters, respectively, where communications only took place after the individual delineation work. The final DBN annotations for the proposed atlas were carried out by applying the consensus of the delineated parcellation maps from the two raters. A voxel was assigned to a region only if both raters labeled the voxel as the same region.
The labeled nuclei include the bilateral Pu, Cau, GPi, GPe, RN, SNr, SNc, NAcc, STN, VP, and Thal. For the nuclei located adjacent to the telencephalon, such as GP (GPi and GPe), Pu, and Cau, the main reference was the T1w atlas, while the QSM atlas was used for auxiliary outline positioning. The Cau and NAcc were mainly delineated in the sagittal view with respect to the division of the lower boundary of cerebrospinal fluid. The lower fiber bundle was considered in Pu parcellation. The Pu, GPi, and GPe were mainly labeled in the axial view and are closely attached. The QSM atlas was mainly referenced in the labeling work of iron-rich DBN, including the STN, RN, VP, and SN (SNr and SNc). The parcellation of STN was mainly defined in the axial view with the spindle shape. The parcellation of spatially close SN was defined with respect to the STN, following the guidance of previous research and histological atlases (Behrens et al. 2003;He et al. 2017;Morel et al. 1997).

Construction of thalamic subnuclei annotations based on high-resolution QSM image and label fusion method
The process of thalamic subnuclei parcellation is shown in Fig. 2. First, a high-resolution (h-reso) QSM image (0.4 × 0.4 × 0.4 mm 3 ) was scanned using a 7T scanner using the parameter presented in "MR Image acquisition and reconstruction". Then, a detailed thalamic subregion parcellation map was created by the two radiologists referring to Schaltenbrand-Wahren's stereotactic atlas, with the assistance of the border-enhancing technique. The h-reso QSM image was then aligned to the low-resolution (l-reso) individual QSM images for compensating the single individual image bias, yielding transformations = {R 1 , R 2 , ⋯ , R N } . For registration detail, l-reso images were first up-sampled to the size of the h-reso image to reduce the influence of resolution difference between scans, and then they were registered to the h-reso image. Then, the thalamic subnucleus parcellation map was correspondingly propagated into each l-reso individual space via applying the inverted deformation fields −1 = {R −1 1 , R −1 2 , ⋯ , R −1 N } and down-sampled to the size of l-reso individual images. This is described in the blue part of Fig. 2. Later, for each individual l-reso QSM image, a detailed thalamic parcellation map was created. Then, based on the deformation fields Φ = { 1 , 2 , ⋯ , N } from QSM atlas construction, joint label fusion (Wang et al. 2012) was performed on each individual thalamic parcellation map. Finally, the 100 individual thalamic subnuclei parcellation Fig. 2 The pipeline of creating the thalamic subnuclei in atlas space. The low-resolution individual QSM images are warped to the highresolution QSM image space, and then the parcellation map is inversely warped to each individual's image, in order to form thalamic parcellation. Finally, the detailed thalamus parcellation map for atlas space is generated by joint fusion maps were fused to form a 64-ROI parcellation of the bilateral thalamus in the 3T atlas space. This is described in the green part of Fig. 2.

Comparison of image contrast inside and outside DBN across different atlases
High image contrast between DBN subregions and the surrounding brain tissues facilitates delineation for DBN subregion boundaries. We performed a quantitative analysis of contrast using the border-enhancing method proposed in "DBN border-enhancement and manual parcellation". We selected several typical image patches containing rich DBN regions and apply border-enhancing operation. The contrast difference can be boosted by the border-enhanced patches. Then, we extracted profile lines across the border-enhanced DBN-rich patches and created intensity profiles of atlases to quantitatively compare the DBN image contrasts of atlases built from different strategies. The profile line was chosen to cross through the most DBN in the patches. More intensive variation along the intensity profile suggests higher contrast in image intensity. In addition, we conducted k-means intensity clustering on selected atlas slices to generate a DBN segmentation map in an unsupervised fashion. A more accurate and detailed clustering result with less small-sized noisy clusters suggests a sharper contrast between the DBN regions and the surrounding tissues.

Analysis of using parcellation maps to delineate DBN in individual images
In practice, manual parcellation is usually needed for individual brain image DBN delineation, especially for smallsize nuclei. Wrapping the predefined parcellation map to the individual images is a more efficient way of DBN delineation. The process includes registering the individual images to the atlas and then applying the inverse warping of the atlas-space parcellation map back to the individual images. We conducted a DBN segmentation test across atlases to assess the accuracy of creating parcellation maps for individual QSM scans. We assumed that the DBN parcellation of our atlas is more accurate than those of the previous atlases.
For evaluation, we computed the Dice similarity coefficient (DSC) and the Hausdorff distance (HD) scores between the inverse warping of the atlas-spaced parcellation maps and the manual delineations of DBN in 6 individuals.
DSC measures the overlap proportion between the warped parcellation and the GT, ranged from 0 to 1. A DSC closer to 1 indicates that the warped parcellation overlaps more with the GT, thus we prefer larger DSCs. The DSC between parcellations A and B is defined as The HD measures the maximal element distance between the warped parcellation and the GT. An HD with smaller value indicates the warped parcellation is generally robust which is bounded more compactly to the GT. The HD between parcellations A and B is defined as Cau, Pu, GPe, GPi, NAcc, SNc, SNr, RN, VP, and STN are the nuclei of interest. The following atlases and parcellations were selected for comparison with ours: Colin-27 atlas (Holmes et al. 1998) with AAL3 parcellation (version 3.1, updated in June 2020) (Rolls et al. 2020), Zhang's longitudinal QSM atlas with its parcellation (Zhang et al. 2018), ICBM 152 2009c Nonlinear Asymmetric T1w atlas (Fonov et al. 2011) with CIT168 subcortical reinforcement learning parcellation (Pauli et al. 2018), AHEAD ICBM-aligned QSM template with its atlas (Alkemade et al. 2020). Based on the different nuclei included and the different imaging properties of the atlases, we categorized the compared atlases into three groups: the coarse-grained DBN parcellation group (the coarse group), the fine-grained DBN parcellation group (the fine group), and the high-resolution DBN parcellation group (the h-reso group).
The coarse group includes Colin-27 atlas with AAL3 parcellation and Zhang's atlas with its proposed parcellation. The Colin-27 T1w and Zhang's QSM atlases were both acquired on 3T scanners with a voxel resolution of 1 × 1 × 1 mm 3 . Colin-27 is in the MNI space, and Zhang's atlas is in the native atlas space. In the coarse group, the nuclei of Cau, Pu, GP, NAcc, and SN were involved for comparison.
The fine group includes ICBM 152 2009c Nonlinear Asymmetric T1w atlas with CIT168 parcellation in comparison. The parcellation of CIT168 is a probability map, and we sampled it to a binary map with a threshold of 0.5, as the original publication suggested. The left-and righthemisphere nuclei are differently marked, for the depiction of the left-and right-hemisphere nuclei are not separated in the original parcellation. The comparison of this group is still intra-modal, respectively, with our QSM atlas and ICBM T1w atlas. All DBN depicted in our work including Cau, Pu, GPe, GPi, NAcc, SNc, SNr, RN, VP, and STN were involved for comparison.
The h-reso group only includes the AHEAD atlas in comparison. The AHEAD atlas is created on multi-modal 7T individual QSM scans with the voxel resolution of 0.5 × 0.5 × 0.5 mm 3 . DBN of GPe, GPi, SN, RN, and STN were involved for comparison.

Application: revealing the PD-related pathology abnormality of iron content in DBN
In this section, we applied the novel-proposed QSM atlas and DBN parcellations to identify the PD-related iron content abnormality, especially in deep brain regions. First, we performed a whole-brain voxel-wise statistical comparison of susceptibility value between PD patients and normal controls (NC). All individual QSM scans were warped to the MNI space by co-registration to the proposed QSM atlas. A smoothing operation was conducted by applying a 3D Gaussian kernel with a standard derivation of 3mm , to further prevent the potential registration error. Then, the statistical examination was performed by a twosample t-test using the randomize utility implemented in FSL (Winkler et al. 2014;Woolrich et al. 2009) with 5000 permutations and threshold-free cluster enhancement. The significance threshold was set at P< 0.05 with family-wise error (FWE) corrected. Before the statistical examination, the covariates of age and gender were first tested whether to be irrelevant across the groups by examining the null hypothesis that the mean values of the tested covariates showed no significant differences between groups. Then, the two groups of scans were tested with the null hypothesis of voxel-wise group medians of susceptibility value were equal. If covariates showed a significant difference between the groups, i.e., the null hypothesis regarding covariates cannot be rejected, the covariates were added as additional variables for a covariate-adjusted test. The constructed DBN parcellation map was applied to identify the PD-related nuclei and thalamic subnuclei that showed significant value in the permutation test.

Results
This section first illustrates the atlas construction performance conducted based on the hybridization strategy. Then, the comparison of image contrast among the atlases is presented. Thirdly, it shows the manual segmentation of the DBN and the thalamic subnuclei on the atlases. The evaluation of individual image DBN parcellation accuracy and application in PD iron content analysis are presented in the last part.

Atlas construction
The hybrid, T1w, and QSM atlases constructed in the purposed hybrid-guided fashion (abbreviated as guided atlases) are shown in Fig. 3, together with the T1w atlas and the QSM atlas built without image hybridization guidance (marked as direct atlases). As illustrated, the T1w atlas shows a clear contrast between gyri and sulci, and the QSM atlas shows clear boundaries of the iron-rich DBN. The T1w-QSM hybrid atlas combines the advantages of the T1w and the QSM atlases and shows high image contrast in both the cortices and the DBN. Compared with the T1w atlas, the depiction of the DBN, such as the thalamus and SN, is greatly improved in the hybrid atlas.

Contrast analysis between modalities and construction strategies of atlases
We compared the intensity contrast among the constructed atlases. The evaluation of image contrast are intensity profile analysis of the border-enhanced display and image segmentation using k-means intensity clustering, as purposed in "Comparison of image contrast inside and outside DBN Fig. 3 Qualitative comparison of the hybrid-guided atlases (guided atlases, ours) and the directly constructed atlases (direct atlases). All atlases are presented in the MNI coordinate space. QSM images are presented using susceptibility of [−0.1, 0.1] ppm across different atlases". We compared the proposed T1w, QSM, multi-modal hybrid-guided atlases, and the ICBM T1w atlas as a reference. Specially, two image patches containing the DBN of interest in the basal ganglia and thalamus are extracted for the intensity profile analysis.
The border-enhanced display and the intensity profile results are shown in Figs. 4 and 5, respectively. The profile line is a specific row on the border-enhanced atlases crossing DBN in the image, and is marked as the black dashed line in Fig. 4. The intensity profile result is the voxel intensity value along the profile line. A large fluctuation in intensity indicates high contrast between the nuclei and the surrounding brain tissues. Fig. 4, both T1w atlases (ICBM and our guided T1w) provide relatively clear boundaries of the nuclei of Cau, Pu, GP, as well as a coarse border of Thal. The division between GPe and GPi can also be observed on both T1w atlases. Our QSM atlas depicts the nuclei that are hard to be observed in T1w atlases, such as the RN and SN. The contrast of borders of other DBN, such as Cau, Pu, and GP (GPe and GPi), are also enhanced. The depiction of the thalamic subnuclei boundary is also clearer on the QSM atlas than T1w atlas. However, the QSM atlas displays poorer contrast in cortical regions than the T1w atlases. The hybrid atlas combines the advantages of the T1w and QSM atlases, displaying a sharp contrast in both cortical and subcortical regions. Both the T1w and the hybrid atlases have similar cortical visual features, and the typical visual features of subcortical ROIs, such as RN and thalamic subnuclei, in the QSM atlas, are also preserved on the hybrid atlas.

Results of intensity profile analysis As shown in
According to the intensity profile result shown in Fig. 5, the same modalities present a similar level of image contrast. The result shows great subcortical sharpness of the QSM atlas, indicated by the largest intensity variations (approximately [0, 252]), and thus retains the clearest features of the DBN anatomical structures. Both T1w atlases show a much smaller variation (approximately [105,150]), although with a similar variation trend as in the QSM result. The hybrid atlas shows a medium intensity variation (approximately [40, 200]). Notably, the result of the hybrid atlas properly combines the result features from both T1w and QSM atlases, showing the same trend as the T1w atlas and retaining the twists and turns of the QSM atlas. Fig. 6, where the generated rough parcellation maps are overlapped with the tissue intensity image. We set k = 4 as the optimal parameter after multiple attempts. The worst segmentation is observed on the guided T1w atlases, where the gray matter and white matter tissues cannot be accurately separated. Conversely, the segmentation on the two QSM atlases accurately delineates the typical DBN, such as Cau and Pu, and thalamus is divided into two parts. However, the white matter tissue and cerebrospinal fluid cannot be accurately segmented on the QSM atlas. The segmentation result of the hybrid atlas outperforms those of the T1w and QSM atlases, where the gray matter is more accurately delineated from the surrounding white matter tissue and cerebrospinal fluid. Furthermore, the delineation of DBN on the hybrid atlas is finer than that on the T1w atlas. Similar segmentation performance compared with those of guided T1w and QSM atlases can be observed on the directed T1w and QSM atlases.

DBN parcellation map
Both the original atlases and the border-enhanced display of the multi-modal atlases (the T1w and QSM atlases) are referred to for the parcellation of the DBN so that we could make full use of the features on the multi-modal atlases. Figure 7 shows multi-modal atlases where delineation of the nuclei of interest are referred. Three typical slices with multi-modal atlases are shown in group, to clarify which image modality is preferred in the manual delineation of different DBN. Meanwhile, manually delineated curves are overlapped on top of the evidencing boundary in nuclei delineation (single hemisphere only for slices in axial view). Note that the curve of each nucleus is delineated on the exactly referred atlas in the nuclei delineation step, and directly translated to the intra-group atlas samples. For example, the Pu is directly delineated on the borderenhanced display of the T1w atlas (Fig. 7a). The GP (GPe and GPi) are delineated on the original QSM atlas since the borderline is not with sharp enough contrast on the T1w or the border-enhanced QSM displays (Fig. 7b, blue and  Fig. 4 are enlarged and plotted to demonstrate image con-trast in the ROI. Image intensities of all Laplacian map patches were scaled to [0,255] in advance red). The Cau and NAcc are depicted on the sagittal view of the original T1w atlas, for only the T1w atlas most clearly shows the NAcc. The Cau is first depicted with reference to the NAcc, as the NAcc is positioned to the inner boundary of the Cau (Fig. 7c). Depiction of other nuclei in the basal ganglia, including the RN, SN, and STN, is performed on the original QSM atlas since the border-enhanced scan is unable to provide more competitive contrast (Fig. 7c). The STN is delineated on the original QSM atlas. The SNr and SNc are depicted on both T1w and QSM atlases since the outer border of the SNr and the inner border of the SNc are inconsistent across modalities. To solve the inconsistency, the SNr was first depicted with reference to the position of the STN, and the SNc is then depicted with reference to the position of the SNr. RN is depicted on the QSM atlas. The references to nuclei depiction are summarized in Table 1. Figure 8 shows the DBN parcellation result, including a 3D rendering view (upper panel) and 2-D slice views covered on the constructed T1w atlas (lower panel). The zoomed-in view of the regions of interest is also shown for closer visualization. The volumes of the nuclei calculated from the parcellation are shown in Table 2.

Accuracy of using parcellation map to delineate DBN in individual images
The comparison results for assessing atlas-guided DBN segmentation accuracy in the individual QSM images are shown in Fig. 9. To make the illustrations of the experiment result clear, we categorized the DBN into two groups, the easy segmentation group and the hard segmentation group. Nuclei with relatively large size and sharp contrast with surroundings are categorized as the easy group. Segmentation in the easy group mainly achieves relatively more accurate (higher mean DSC and lower mean HD) and more stable (small variance) performance across all the atlases in comparison. Contrarily, nuclei with relatively small size and blurred contrast with surroundings are categorized into the hard group. Segmentation performance of DBN in the hard group mainly turns out to be less accurate and more likely to appear as outliers or large variance in evaluation metrics. In the DBN segmentation accuracy evaluation, we consider the overall accuracy superior to the robustness. Relatively unstable results with better overall accuracy are considered better than stable but inaccurate results. The category of each nucleus is considered both from visualization and performance perspectives. In detail, we categorize the nuclei of Pu, GP (GPe), and RN as the easy segmentation group, while the nuclei of NAcc, Cau, STN, SN (SNr/SNc), GPi, and VP are categorized into the hard segmentation group.
Our atlas outperforms the previous atlases by 5% to 10% in both evaluation metrics (DSC and HD). The means of DSCs using our parcellations are higher than 0.8, ranging from 0.73 to 0.90. In contrast, the DSCs using the previous atlases are mostly lower than 0.8, ranging from 0.6 to 0.7. The lower HD values also indicate better performance on DBN delineation using our parcellation than the previous parcellations.
For the coarse group, our delineation of all the nuclei achieves higher DSCs compared with that of the previous atlases. The DSCs are high with small variance for the GP,  Pu, Cau, and RN. The DSCs of SN and NAcc are also higher than the other atlases. Our atlas also shows a relatively shorter HD than Zhang's and AAL3 atlases, with lower mean value and variance (e.g., the distance of our results versus that of Zhang's atlas). Both of the other atlases show suboptimal results for bilateral Cau, while our atlas achieves a low mean in the distance with slightly high variance, compared with significantly higher means in the distance of the other atlases. Notice that NAcc is not included in Zhang's atlas, and thus the corresponding results of both metrics are presented as 0.
For the comparison of the h-reso group, the delineation using our parcellation outperforms the AHEAD atlas from both DSC and HD perspectives with better mean value and lower variance. For difficult nuclei such as STN and SN, our atlas still shows competitive performance with both reasonable value and variance not out of the ordinary. Figure 10 shows the axial and coronal views of the manually annotated thalamic subnucleus parcellation map overlays on a high-resolution (7T) QSM image, with the histology reference shown on the right side. The QSM image with a voxel resolution of 0.4×0.4×0.4 mm 3 provides clear contrast among thalamic subnuclei for annotation. Table 3 (shown in the appendix) presents the list of the total 32 thalamic subnuclei on each side and the corresponding color codes. Figure 11 illustrates two representative slices of the 3T QSM atlas with fused parcellation map overlays on the top, with the histology reference shown on the right side. The fused MNI space subnucleus labels of the thalamus show more smooth boundaries than that in the individual 7T guidance image. The label fusion strategy facilitates histology-consistent thalamus subnucleus annotation in the atlas space. In the final parcellation map, we can precisely identify crucial thalamic subnucleus such as the Ventrointermedius (V.i.m., highly related to PD), Dorsomedialis (M, related to schizophrenics), Anteroprincipalis (A.pr, a vital component of the hippocampal system for episodic memory).

DBN with PD-related abnormality in iron content identified on the atlas
The map revealing the group differences in brain tissue susceptibility between the PD patients and the normal controls indicated by the voxel-wise unpaired t-test is shown in Fig. 12 using ITK-SNAP (Yushkevich et al. 2006). The voxels with FWE-corrected P-values less than 0.05 are overlaid on the MNI space QSM atlas. For focusing attention on DBN regions, we masked the non-DBN susceptibility variations out of our significance map. The boundaries of DBN are further overlaid on the significance map to identify the DBN with disease-related alternations in iron content. Figure 13 shows a representative slice with detailed DBN parcellation. The abnormalities in iron content are found in the bilateral Cau, Pu, GP, Thal, and RN, as well as the right STN nucleus and SN (both SNr and SNc). The volume of regions with abnormality in QSM in the left hemisphere is slightly larger than that in the right hemisphere.
The abnormality in the thalamus mainly distributes from the ventral region to the central medial region, approximately from the coordinate of −5 mm to 13mm along the Z-axis, shaped like a spindle with its most part on the ventral   Our results are consistent with the previous reports of the PD-related alterations in iron content. For example, a study using the voxel-wise permutation test revealed QSM abnormality in the SN and STN . The abnormality in SN is consistently found in PD patients as revealed by a systematic review that reported 30 out of 33 (90.9%) studies indicated significant PD-related QSM abnormality in the SN (Ravanfar et al. 2021). This review also suggested a wide range of brain areas, including the Cau, Pu, and Thal, with higher susceptibility in PD patients with cognitive impairments.

Discussion
In the present study, we first introduce the MuSus-100 hybrid image-guided MNI space QSM atlas construction and the DBN parcellation based on these atlases. Quantitative analysis indicates that the T1w-QSM hybrid atlas better depicts the cortical and subcortical structures than the T1w atlas or QSM atlas alone. Furthermore, the MuSus-100 atlas outperforms the current existing atlases in the accuracy and fineness of DBN segmentation in individual images. In addition, a fine parcellation of the thalamic subnuclei is introduced via a label fusion strategy to our 3T atlas using a manual parcellation on a 7T QSM image as guidance.

Quality of constructed atlases and parcellation maps
For the quality of the native space atlas, the construction steps including image registration and projection processes are well controlled by introducing suitable processing steps, including the hybrid-image-guided group-wise registration, the average template sharpening update, and the T1w atlasguided MNI space projection. In particular, the hybrid image combines the advantages of both the T1w image and the QSM image, which results in a more accurate registration in both the cortical and subcortical structures across the samples. On the other hand, our atlases are constructed from 100 young adult volunteers with ages ranging from 19 to 29. The age distribution is approximately identical to that of typical MNI space atlases of Colin-27 (28-year individual) (Holmes et al. 1998) and ICBM 152 atlas series (18.5-43.5) (Mazziotta et al. 2001;Fonov et al. 2009Fonov et al. , 2011. With a  Significance map ( P FWE < 0.05 ) revealing the difference in tissue susceptibility between PD patients and healthy controls indicated by a voxel-wise permutation test with age-correction. The between-group comparison heat map result is overlaid on the constructed MNI space QSM atlas NAcc are relatively hard to depict due to the coarse boundary and fuzzy local features.
For the difficult DBN, the performance of both our parcellation and those based on the atlases for comparison in delineating DBN in the individuals are not ideal. The propagation of these atlases results in relatively low mean DSCs. Moreover, the outliers in the result of HD are more likely to occur for these nuclei. Nevertheless, our atlas still outperforms those for comparison, with fewer outliers and equivalent DSC and HD scores for the nuclei difficult to depict. Thus, our atlas still provides relatively satisfying results in the nuclei with fuzzy surrounding information.
For the segmentation result of the fine group, the CIT168 atlas has only structural imaging atlases, T1w and T2w. The susceptibility differences between DBN and surroundings provide intrinsic advantages of QSM-based atlas in the DBN segmentation. The reason for choosing the CIT168 atlas for comparison is that only the CIT168 atlas can offer this fine-grained level DBN parcellation in the MNI space. We would like to compare the projection accuracy of each nucleus defined in the MuSus-100 atlas. The result shows that our atlas achieved state-of-the-art performance at this fine-grained level of DBN parcellation.
The small size of the validation image set makes drawbacks to the atlas evaluation experiments. Since it is timeconsuming to generate accurate DBN labels, we only use 6 labeled individual images for testing in the experiment. Therefore, even only one imperfect segmentation result will cause a severe effect on average value and outlier or large variance in the result plots. For example, for one test subject, the HDs between our segmentation and GT in the CauR is larger than that in other subjects (about 40 voxels compared with 15 voxels on average). This is generally because Cau is the largest one among all the nuclei for segmentation tests, and its long tail may bring difficulty for accurate registration. Thus, a minor registration error (which is hardly indicated by DSC) can cause a large HD value. Fig. 13 A representative slice of QSM atlas with significance map ( P FWE < 0.05 ) and DBN boundaries overlaid on top. The nuclei with extensive coverage of significant voxels are indicated by arrows. The names of thalamic subregions are shown in the sheet displayed on the right similar age distribution, the imaging result of brain structure and tissue is similar between our T1w atlas and the reference MNI space T1w atlases. Thus, accurate projection of the native spaced atlases to the MNI space is feasible, resulting a accurate projection result for our atlas in the MNI space.
Both visual inspection and contrast analysis demonstrate high quality for the constructed atlases in delineating the DBN of interest. By direct visual inspection, several smallsize DBN and the DBN anatomical subregions can be clearly identified, which are usually invisible or not clear in other atlases. For instance, the division of GPe and GPi can be observed even in our T1w atlas. We performed the k-means intensity clustering on selected image patches, and the border-enhanced displays to assess the image quality for constructed atlas. Both analyses show that the QSM atlas provides a sharp contrast in DBN regions, which improves the delineation performance for nuclei of interest. In addition, both the constructed T1w and QSM atlases and the borderenhanced display are considered references to maximally improve the accuracy of DBN parcellation.

Accuracy of the atlas-based DBN parcellation of individual images
The results show that our atlas outperforms those of the existing DBN parcellations based on the other T1w and QSM atlases both in the accuracy and stability of the work of propagating the parcellation map to the individual QSM images. For all nuclei, the delineation generated by our parcellation outperforms those for comparison in both mean and variance, which indicates that our result provides higher accuracy and stability. Since the difficulty of the depiction of nuclei may vary due to different local neuroimaging contrast and inter-individual variations, the nuclei easier to depict may achieve higher accuracy and stability. In our present study, the GP, RN, Cau, and VP are relatively easy to depict for they have sharp boundaries, while the STN, SN, and

Application prospects of the atlases and parcellation maps
Accurate and fine parcellations of subcortical nuclei, particularly the nuclei within the basal ganglia and thalamus, are greatly needed in the studies of neurological and psychiatric. Specific types of symptoms of neurological and psychiatric diseases may be associated with abnormal structural or functional features in specific DBN. Furthermore, deep brain stimulation (DBS) targeting the DBN has been used to treat neurological and neuropsychiatric disorders in patients who could not receive sufficient effect from non-invasive therapies (Williams and Okun 2013). Studies examining the relationship between the anatomical location of DBS and the treatment effect have been carried out recently to identify the optimal location of stimulation (Bot et al. 2018;Matsumoto et al. 2016). Thus, the faithful depiction of DBN based on neuroimaging is crucial for the studies of disease pathology and target refining for deep brain stimulation.
Our atlases derived from healthy young volunteers (around 25 years old) may serve as a representation of a healthy young population. It is known that human brain goes through dynamic age-dependent variations in brain structure and tissue composition (Fjell et al. 2013;Giorgio et al. 2010;Cabeza et al. 2018). We thus expect that our atlas can provide the best potential applications in studying neurological disorders within the age range 18-40 years old. Compared with elder population involved in the longitudinal studies on QSM alterations (Zucca et al. 2006;Ward et al. 2019;Peterson et al. 2019), our atlas would not achieve full representativeness. For example, Zhang et al. (2018) performed longitudinal analysis on the quantitative susceptibility changes using fine labeled QSM and T1w scans. For individual groups with different ages (e.g., the elder), gender (e.g., female-dominated), and diseases (e.g., PD), atlases built specifically on those groups will provide a better representation.

Mapping of thalamus subnuclei and its application prospects
As reported by a series of recent studies, QSM based on the GRE sequence is superior in delineating the thalamic subnuclei (Deistung et al. 2013;He et al. 2020;Li et al. 2019a). Here, we for the first time provide a fine parcellation of the thalamic subnuclei based on QSM from a 7T MRI, which identifies 32 subnuclei of the thalamus on each side confirmed with the Schaltenbrand and Wahren histology atlas. Due to the high expense of manual segmentation, we manually delineated the thalamic subnuclei on only one 7T image as the GT. This might introduce individual bias. To solve this issue, we projected the parcellation onto all the 3T individual images by nonlinear registration and afterward fused the individual parcellations in the atlas space to achieve the final parcellation of thalamic subnuclei. Such a process should be able to guarantee our parcellation unbiased.
Using the fine parcellation map of the thalamus, we can precisely identify critical thalamic subnuclei such as the Ventrointermedius (VIM, highly related to tremor), Dorsomedialis (M, related to schizophrenics), and Anteroprincipalis (A.pr, a vital component of the hippocampal system for episodic memory). Thus, the proposed thalamic parcellation map in standard space can become a powerful tool for understanding related neurological disease progression and providing insights for the treatments of diseases related to abnormalities. This fine atlas of thalamic subnuclei could also provide a visual reference to confirm the accuracy of indirect targeting, which may potentially help to improve the treatment effect of the surgeries on thalamic subnuclei. Meanwhile, in neurosurgery experiences, the proposed thalamic parcellation map allows positioning for critical targets in DBS surgeries using T1w and T2w scans. In addition, the fine parcellation also boosts research on subthalamic connectomics, which offers fine-targeted ROIs for finding the relationship between related diseases and the abnormal ROI connections in a more precious view in data-abundant 3T space.

Conclusions
In the present study, we introduce an MNI space human brain QSM atlas and parcellation map, named MuSus-100 atlas. The QSM atlas is derived from the T1w, QSM, and GRE magnitude images of 100 healthy young volunteers via group-wise registration guided by the T1w-QSM linear hybridization of the individual images. The hybrid image highlights the representation of subcortical nuclei, and meanwhile, preserves the advantages of clear anatomical tissue structure in T1w images. The constructed MNI space QSM atlas provides sharp contrast for facilitating the delineation of the DBN of interest. Ten highly interested DBN are depicted based on atlases of multiple modalities. In addition, a fine thalamus parcellation map based on the 7T individual QSM is introduced into MNI space, providing a detailed histology-consistent thalamic subnucleus parcellation for MNI space images acquired from a 3T scanner. The constructed parcellation maps clearly delineate DBN of interest in the standard MNI space, which can help the ROI-based quantitative analysis of iron-rich subcortical nuclei in future studies. The proposed atlas has the potential value to facilitate the research in the field of neurological and psychiatric disease progress and DBS surgical target refining.