Multicontrast anatomical subcortical structures parcellation
Abstract
The human subcortex is comprised of more than 450 individual nuclei which lie deep in the brain. Due to their small size and close proximity, up until now only 7% have been depicted in standard MRI atlases. Thus, the human subcortex can largely be considered as terra incognita. Here, we present a new opensource parcellation algorithm to automatically map the subcortex. The new algorithm has been tested on 17 prominent subcortical structures based on a large quantitative MRI dataset at 7 Tesla. It has been carefully validated against expert human raters and previous methods, and can easily be extended to other subcortical structures and applied to any quantitative MRI dataset. In sum, we hope this novel parcellation algorithm will facilitate functional and structural neuroimaging research into small subcortical nuclei and help to chart terra incognita.
Introduction
Subcortical brain structures are often neglected in neuroimaging studies due to their small size, limited interregional contrast, and weak signaltonoise ratio in functional imaging (Forstmann et al., 2016; JohansenBerg, 2013). Yet, these small and diverse structures are prominent nodes in functional networks (Marquand et al., 2017; Ji et al., 2019), and they undergo pathological alterations already at early stages of neurodegenerative diseases (Andersen et al., 2014; Koshiyama et al., 2018). Deep brain stimulation surgery, originally performed to reduce motor symptoms in essential tremors, is now a promising therapeutic option in later stages of Parkinson’s disease and movement disorders, as well as refractory psychiatric illnesses in obsessivecompulsive disorder, anorexia, or depression (Forstmann et al., 2017; Mosley et al., 2018). Evolutionary genetics even uncovered that in modern humans, Neanderthalinherited alleles were preferentially downregulated in subcortical and cerebellar regions compared to other brain regions (McCoy et al., 2017), suggesting these structures to be essential in making us specifically human.
Despite their importance, these areas are particularly difficult to image. Furthermore, the size, shape, and location of these brain regions changes with development and aging (Fjell et al., 2013; Keuken et al., 2013; Yeatman et al., 2014; Herting et al., 2018). Experiencebased plasticity continuously remodels myelin (Tardif et al., 2016; Hill et al., 2018; Turner, 2019), iron and other magnetic substances accumulate with age or pathology (Andersen et al., 2014; Zhang et al., 2018), both bringing changes in the MRI appearance of subcortical regions with diverse tissue characteristics (Draganski et al., 2011; Keuken et al., 2017).
Thus, mapping the structure and function of the subcortex is a major endeavor as well as a major challenge for human neuroscience. Extensive work available from animal brain models unfortunately does not translate in a straightforward way to human subcortical anatomy nor does it shed much light on its involvement in human cognition (Steiner and Tseng, 2017). Besides serious difficulties in obtaining adequate measures of subcortical neural activity in functional MRI (de Hollander et al., 2017; Miletić et al., 2020), atlases and techniques for labeling accurately and reliably individual subcortical structures have also been scarce (Frazier et al., 2005; Chakravarty et al., 2006; Ahsan et al., 2007; Yelnik et al., 2007; Qiu et al., 2010; Patenaude et al., 2011), typically labeling the thalamus, striatum (or its subdivision into caudate and putamen), and globus pallidus (internal and external segments combined), sometimes the amygdala. However, recent advances in anatomical MRI, combining multiple contrasts and/or quantitative MRI mapping and utilizing the higher resolution achievable with 7 Tesla (7T) and above have started to reduce the gap, each mapping a few additional structures or substructures, primarily the ironrich substantia nigra, red nucleus and subthalamic nucleus (Keuken et al., 2013; Xiao et al., 2015; Visser et al., 2016a; Visser et al., 2016b; Wang et al., 2016; Makowski et al., 2018; Ewert et al., 2018; Iglesias et al., 2018; Pauli et al., 2018; Sitek et al., 2019). While these efforts generated valuable atlases, they do not yet enable to identify many subcortical structures in individual subjects. Manual delineation, on the other hand, requires extensive manual labor from highly trained experts which cannot be easily applied to large cohorts or clinical settings.
Here, we propose a new automated parcellation technique to identify and label 17 individual subcortical structures of varying size and composition in individual subjects, based on a large quantitative 7T MRI database (Alkemade et al., 2020), using quantitative maps of relaxation rates R1 and R2* (1/T1 and 1/T2*, respectively) and quantitative susceptibility maps (QSM) as anatomical contrasts. The algorithm, named Multicontrast Anatomical Subcortical Structure Parcellation (MASSP), follows a Bayesian multiobject approach similar in essence to previous efforts (Fischl et al., 2002; Eugenio Iglesias et al., 2013; Visser et al., 2016a; Garzón et al., 2018), combining shape priors, intensity distribution models, spatial relationships, and global constraints. The main innovation of our approach is to explicitly estimate interfaces between subcortical structures based on a joint model derived from signed distance functions. Modeling interfaces in addition to the structure itself provides a rich basis to encode relationships and anatomical knowledge in shape and intensity priors. A voxelwise Markovian diffusion regularizes the combined priors for each defined interface, lowering the imaging noise. Finally, the voxelwise posteriors for the different structures and interfaces are further combined into global anatomical parcels by topology correction and region growing taking into account volumetric priors, which regularizes parcellation results further in smaller nuclei with low or heterogeneous contrast. To validate the results from this new method, in a thorough comparison with expert manual labeling, we show that the proposed method provides results very close from manual raters in many structures and exhibit reasonable bias across the adult lifespan. The method can easily be extended to new structures, can be applied to any quantitative MRI dataset and is available in Open Source as part of Nighres (Huntenburg et al., 2018), a neuroimage analysis package aimed at highresolution neuroimaging.
Results
The MASSP parcellation method presented here has been trained to parcellate the following 17 structures: striatum (Str), thalamus (Tha), lateral, 3rd and 4th ventricles (LV, 3V, 4V), amygdala (Amg), globus pallidus internal segment (GPi) and external segment (GPe), SN, STN, red nucleus (RN), ventral tegmental area (VTA), fornix (fx), internal capsule (ic), periaqueductal gray (PAG), pedunculopontine nucleus (PPN), and claustrum (Cl), see Figure 1. These structures include the most commonly defined subcortical regions (Str, Tha, Amg, LV), the main ironrich nuclei (GPi, GPe, RN, SN, STN), as well as smaller, less studied areas (VTA, PAG, PPN, Cl), white matter structures (ic, fx), and the central ventricles (3V, 4V).
MASSP uses a data set of ten expert delineations as a basis for its modeling. From the delineations, an atlas of interfaces between structures, shape skeletons, and interface intensity histograms are generated, and used as prior in a multiplestep noniterative Bayesian algorithm, see Figure 2 and Materials and methods.
Validation against manual delineations
In a leaveoneout validation study comparing performance with the manual delineations, MASSP performed above 95% of the level of quality of the raters for Str, Tha, 4V, GPe, SN, RN, VTA, ic in terms of Dice overlap, the most stringent of the quality measures (see Figures 3 and 4 and Table 1). Several of the smaller structures have lower overlap ratios likely due to their smaller size (GPi, STN, PAG, PPN). Structures with an elongated shape (fx, Cl) remain challenging, due to the fact that small differences in location can substantially reduce overlap (Bazin et al., 2016). Despite these challenges, when comparing the dilated Dice scores, all structures were above 75% of overlap, with most reaching over 90% of the manual raters ability. Note that the Dice coefficient is very sensitive to size, as smaller structures will have lower overlap ratios for the same number of misclassified voxels. The dilated Dice coefficient is more representative of the variability regardless of size, as the smaller structures can reach high levels of overlap, both in manual and automated parcellations (see Table 1). The average surface distance confirms these results, showing values generally between one and two voxels of distance at a resolution of 0.7 mm, except in the cases of Amg, LV, fx, PPN, and Cl. These structures are generally more variable (LV), elongated (fx, Cl), or have a particularly low contrast with neighboring regions (Amg, PPN).
Comparison to other automated methods
To provide a basis for comparison, we applied other freely available methods for subcortical structure parcellation to the same 10 subjects. MASSP performs similarly to or better than Freesurfer, FSL FIRST and a multiatlas registration using ANTs (see Table 2). Multiatlas registration provides high accuracy in most structures as well, but is biased toward underestimating the size of smaller and elongated structures where overlap is systematically reduced across the individual atlas subjects. Multiatlas registration is also quite computationally intensive when using multiple contrasts at high resolution. Finally, MASSP provides many more structures than Freesurfer and FSL FIRST, and can be easily applied to new structures based on additional manual delineations.
Application to new MRI contrasts
Quantitative MRI has only become recently applicable in larger studies, thanks in part to the development of integrated multiparameter sequences (Weiskopf et al., 2013; Caan et al., 2019). Many data sets, including largescale open databases, use more common T1 and T2weighted MRI. In order to test the applicability of MASSP to such contrasts, we obtained the testretest subset of the Human Connectome Project (HCP, Van Essen et al., 2013) and applied MASSP to the 45 preprocessed and skullstripped T1 and T2weighted images from each of the two test and retest sessions. While performing manual delineations on the new contrasts would be preferable, the model is already rich enough to provide stable parcellations. Testretest reproducibility is similarly high for MASSP and Freesurfer, and are generally in agreement, see Figure 5 and Table 3.
Biases due to atlas size
A common concern of brain parcellation methods is the risk of biases, as they are typically built from a small number of manual delineations. Our data set is part of a large scale study of the subcortex, for which we obtained manual delineations of the STN, SN, RN, GPe, and GPi on 105 subjects over the adult lifespan (18–80 year old, see Alkemade et al., 2020 for details). First, we investigated the impact of atlas size. We randomly assigned half of the subjects from each decade to two groups, and built atlas priors from subsets of 3, 5, 8, 10, 12, 15, and 18 subjects from the first group. The subjects used in the atlas were taken randomly from each decade (1830, 3140, 4150,5160, 6170, 7180), so as to maximize the age range represented in each atlas. Atlases of increasing size were constructed by adding subjects to previous atlases, so that atlases of increasing complexity include all subjects from simpler atlases. Results applying these atlases to parcellate the second group are given in Figure 6. As in previous studies (Eugenio Iglesias et al., 2013; Bazin and Pham, 2008), performance quickly stabilized with atlases of more than five subjects (no significant difference in Welch’s ttests between using 18 subjects or any subset of 8 or more for all structures and measures).
Biases due to age differences
To more specifically test the influence of age on parcellation accuracy, we defined again six age groups by decade and randomly selected 10 subjects from each group. Each set of subjects was used as priors for the five structures above, and applied to the other age groups. Results are summarized in Figure 7. Examining this age bias, we can see a decrease in performance when parcellating subjects in the range of 60 to 80 years of age. The choice of priors seem to have a limited impact, which varies across structures. In particular, using priors from a similar age group is not always beneficial.
Bias on individual measures
Finally, we investigated the impact of this decrease in performance in the estimation of anatomical quantities, see Figure 8. The bias did affect the morphometric measures of structure volume and thickness, but the effects on the local measure of thickness was reduced compared to the global measure of volume. Quantitative MRI averages were very stable even when age biases are present in the parcellations.
For reference, we report structure volumes, thickness, R1, R2* and QSM values estimated from the entire AHEAD cohort for different age groups, extending our previous work based on manual delineations on a different data set (Keuken et al., 2017; Forstmann et al., 2014). Results are given in Table 4, describing average volumes, thickness, and quantitative MRI parameters for young, middleaged, and older subjects for the 17 subcortical structures.
Discussion
Our goal with the MASSP algorithm was to provide a fully automated method to delineate as many subcortical structures as possible on highresolution structural MRI now available on 7T scanners. We modeled 17 distinct structures, taking into account location, shape, volume, and quantitative MRI contrasts to provide individual subject parcellations. Based on our results, we can be confident that the automated parcellation technique performs comparably to human experts, providing delineations within one or two voxels of the structure boundaries (dilated Dice overlap over 75% for all structures, including in aging groups). Results were nearly indistinguishable from expert delineations for eight major structures (Str, Tha, 4V, GPe, SN, RN, VTA, ic), and smaller structures retain high levels of overlap, comparable to trained human raters. This parcellation includes the most commonly defined structures (Str, Tha, SN, RN, STN) with overlap scores comparable to those previously reported (Garzón et al., 2018; Visser et al., 2016a; Eugenio Iglesias et al., 2013; Chakravarty et al., 2013; Patenaude et al., 2011). More importantly, it also includes structures seldom or never before considered in MRI atlases and parcellation methods, such as GPe, GPi, VTA, 3V, 4V, ic, fx, PAG, PPN, Cl. The technique handles structures of varying sizes well, as indicated by dilated overlap and boundary distance. Additional structures can be added, if they can be reliably delineated by expert raters on singlesubject MRI at achievable resolutions. Some enhancement techniques such as building a multisubject template (Pauli et al., 2018) or adding a denoising step (Bazin et al., 2019) may be beneficial. Coregistration to a highprecision atlas as in Ewert et al., 2018 may also improve the initial alignment over the MASSP group average template.
Age biases are present both in expert manual delineations and automated parcellation techniques. Age trajectories in volume and quantitative MR parameters indicate systematic shifts in contrast intensities and an increasing variability with age, associated with changing myelination, iron deposition, and brain atrophy (Draganski et al., 2011; Daugherty and Raz, 2013; Fjell et al., 2013; Keuken et al., 2017). These changes seem only to impact the parcellation accuracy for age groups beyond 60 years and agematched priors did not provide specific improvements, thus indicating that an explicit modeling of age effects may be required to further improve parcellation quality in elderly populations. These results also point to exercising caution when applying automated parcellation methods to study morphometry in elderly or diseased populations, where measured differences may include biases. They also point out that while global volume and local thickness are indeed affected by such biases, quantitative MRI measures are much more robust. Note that this bias is likely present is many automated methods, although they have not been systematically investigated due to the extensive manual labor required. Interestingly, biases also exist in expert delineations: when the size or shape of a structure is refined in neuroanatomical studies, experts may become more or less conservative in their delineations. Automated methods provide a more objective measure in such case, as the source of their bias is explicitly encoded in the atlas prior delineations and computational model. Important applications of subcortical parcellation also include deepbrain stimulation surgery (Ewert et al., 2018), where the number of structures parcellated by MASSP can help neurosurgeons orient themselves more easily, although precise targeting will still require manual refinements, especially in neurodegenerative diseases.
We observed that dilated overlap, that is, the overlap of structures up to one voxel, provided a measure of accuracy largely independent of size, for automated or manual delineations. Imprecision in the range of one voxel in the boundary is to be expected due to partial voluming which impacts Dice overlap. The dilated overlap measure is a better representative of performance and indicates that conservative or inclusive versions of the subcortical regions can be obtained by eroding or dilating the estimated boundary by a single voxel. Such masks may be useful when separating functional MRI signals between neighboring nuclei or when locating smaller features inside a structure. Additionally, the Bayesian estimation framework provides voxelwise probability values, which can also be used to further weight the contribution of each voxel within a region in subsequent analyses.
In summary, our method provides fast and accurate parcellation for subcortical structures of varying size, taking advantage of the high resolution offered by 7T and the specificity of quantitative MRI. The algorithm is based on an explicit model of structures given in a Bayesian framework and is free of tuning parameters. Given a different set of regions of interest or different populations, new priors can be automatically generated and used as the basis for the algorithm. If more MRI contrasts are available, the method can also be augmented to take them into account. The main requirement for the technique is a set of manual delineations of all the structures of interest in a small group of representative subjects. Performance may further improve with the number of included structures, as the number of distinct interfaces increases, refining in particular the intensity priors. In future works, we plan to include more structures or substructures and model the effects of age on the priors. We hope that the method, available in open source, will help neuroscience researchers to include more subcortical regions in their structural and functional imaging studies.
Materials and methods
Data acquisition
Request a detailed protocolOur parcellation method has been developed for the MP2RAGEME sequence (Caan et al., 2019). Briefly, the MP2RAGEME consists of two interleaved MPRAGEs with different inversions and four echoes in the second inversion. Based on these images, one can estimate quantitative MR parameters of R1, R2* and QSM. In this work, we used the following sequence parameters: inversion times TI1,2 = 670 ms, 3675.4 ms; echo times TE1 = 3 ms, TE2,1–4 = 3, 11.5, 19, 28.5 ms; flip angles FA1,2 = 4°, 4°; TRGRE1,2 = 6.2 ms, 31 ms; bandwidth = 404.9 MHz; TRMP2RAGE = 6778 ms; SENSE acceleration factor = 2; FOV = 205×205 x 164 mm; acquired voxel size = 0.70×0.7 x 0.7 mm; acquisition matrix was 292 × 290; reconstructed voxel size = 0.64×0.64 x 0.7 mm; turbo factor (TFE) = 150 resulting in 176 shots; total acquisition time = 19.53 min.
T1maps were computed using a lookup table (Marques et al., 2010). T2*maps were computed by leastsquares fitting of the exponential signal decay over the multiecho images of the second inversion. R1 and R2* maps were obtained as the inverse of T1 and T2*. For QSM, phase maps were preprocessed using iHARPERELLA (integrated phase unwrapping and background phase removal using the Laplacian) of which the QSM images were computed using LSQR (Li et al., 2014). Skull information was removed through creation of a binary mask using FSL’s brain extraction tool on the reconstructed uniform T1weighted image and then applied to the quantitative contrasts (Smith, 2002). As all images were acquired as part of a single sequence, no coregistration of the quantitative maps was required (see Figure 9).
Anatomical structure delineations
Request a detailed protocolManual delineations of subcortical structures were performed by two raters trained by an expert anatomist, according to protocols optimized to use the better contrast or combination of contrasts for each structure and to ensure a consistent approach across raters. The following 17 structures were defined on a group of 10 subjects (average age 24.4, eight female): striatum (Str), thalamus (Tha), lateral, 3rd and 4th ventricles (LV, 3V, 4V), amygdala (Amg), globus pallidus internal segment (GPi) and external segment (GPe), SN, STN, red nucleus (RN), ventral tegmental area (VTA), fornix (fx), internal capsule (ic), periaqueductal gray (PAG), pedunculopontine nucleus (PPN), and claustrum (Cl). Separate masks for left and right hemisphere were delineated except for 3V, 4V, and fx. In the following the algorithm treats each side separately, resulting in a total of 31 distinct structures (see Figure 1).
Anatomical interface priors
Request a detailed protocolIn order to inform the algorithm, we built a series of priors derived from the manual delineations. Each subject was first coregistered to a MP2RAGEME anatomical template built from 105 subjects coaligned with the MNI2009b atlas (Fonov et al., 2011) with the SyN algorithm of ANTs (Avants et al., 2008) using successively rigid, affine, and nonlinear transformations, high levels of regularization as recommended for the subcortex (Ewert et al., 2019) and mutual information as cost function.
The first computed prior is a prior of anatomical interfaces, recording the most likely location of boundaries between the different structures, defined as follows. Given two delineated structures $i,j$, let ${\phi}_{i},{\phi}_{j}$ be the signed distance functions to their respective boundary, that is, ${\phi}_{i}(x)$ is the Euclidean distance of any given voxel to the boundary of i, with a negative sign inside the structure. Then we define the interface ${B}_{ij}$ with the distance function ${d}_{ij}$:
where $\delta $ is a scale parameter for the thickness of the interface. These interfaces functions are not symmetrical, as the intensity inside i next to j is generally different from the intensity inside j next to i. Based on this definition, the prior for a given interface based on N manual delineations is given by:
These probability functions are calculated for all possible configurations including $ii$, which represent the inside of each structure. We thus have a total of ${N}^{2}$ functions, but only a few are nonzero at a given voxel x, and we may keep only the 16 largest values to account for any number of interfaces in 3D (Bazin et al., 2007). Finally, we need to scale the prior to be globally consistent with the priors below by assuming that the 95th percentile of the highest kept $P(x\in {B}_{ij})$ values have a probability of 0.95. The scale parameter $\delta $ is set to one voxel, representing the expected amount of partial voluming. The resulting interface prior is shown in Figure 10A.
Anatomical skeleton priors
Request a detailed protocolNext, we defined priors for the skeleton of each structure, representing their essential shape regardless of exact boundaries (Blum, 1973). As we are mostly interested in the most likely components of the skeleton or medial axis ${S}_{i}$, we follow a simple method to estimate its location:
We define as ${s}_{i}(x)$ the signed distance function of this discrete skeleton, and define prior probabilities as above:
The skeletons are defined inside each structure, which implies $P(x\in {S}_{i})\le P(x\in {B}_{ii})$. To respect this relationship, we scale $P(x\in {S}_{i})$ with the same factor as $P(x\in {B}_{ii})$ but use $\sqrt{P(x\in {S}_{i})}$ when combining probabilities during the estimation stage. The obtained anatomical skeleton priors are given on Figure 10B.
Interface intensity priors
Request a detailed protocolWhile anatomical priors already provide rich information, they are largely independent of the underlying MRI. From the coaligned quantitative MRI maps and manual delineations, we defined intensity priors for every interface $ij$, in the form of intensity histograms to ensure a flexible representation of intensity distributions. Given a quantitative contrast ${R}_{n}(x)$, we built a histogram ${H}_{ij,n}$ for each subject n and interface $ij$. Histograms have 200 bins covering the entire intensity range within a radius of 10 mm from any of the delineated structures. To obtain an average histogram, we combine each histogram with a weighting function ${w}_{n}(x)$ giving the likelihood of the subject’s intensity measurement compared to the group:
where ${\mu}_{R}(x)$ is the median of the ${R}_{n}(x)$ values at x, and ${\sigma}_{R}(x)$ is 1.349 times the interquartile range of ${R}_{n}(x)$. These are robust estimators of the mean and standard deviation, used here to avoid biases by intensity outliers. To further combine the R1, R2*, and QSM contrasts we take the geometric mean of the histogram probabilities: ${H}_{ij}(x)={\prod}_{R}{H}_{ij}{(R(x))}^{1/3}$.
Global volume priors
Request a detailed protocolThe last type of priors extracted from manual delineations are volume priors for each of the structure. Here, we assume a lognormal distribution for the volumes ${V}_{i}$ and simply estimate the mean ${\mu}_{V,i}$ and standard deviation ${\sigma}_{V,i}$ of $\mathrm{log}{V}_{i,n}$ over the subjects.
Voxelwise posterior probabilities
Request a detailed protocolWhen parcellating a new subject, we first coregister its R1, R2*, and QSM maps jointly to the template and use the inverse transformation to deform the anatomical priors into subject space. Then we derive voxelwise posteriors as follows:
Once again we should compute all possible combinations, but due to the multiplication of the priors we can restrict ourselves to the 16 highest probabilities previously estimated. To balance the contribution of the anatomical priors and the intensity histograms, we also need to normalize the intensity priors sampled on the subject’s intensities. We use the same approach, namely assuming that the 95th percentile of the highest kept ${H}_{ij}(x)$ values have a probability of 0.95, separately for each contrast. The voxelwise parcellation and posteriors obtained are shown in Figure 11A.
Markovian diffusion
Request a detailed protocolThe voxelwise posteriors are independent from each other and do not reflect the continuous nature of the structures. The next step is to combine information from neighboring voxels. We define a sparse Markov Random Field model for the posteriors:
with $P(y\sim xR)={\prod}_{R}\mathrm{exp}{(R(y)R(x))}^{2}/2{\sigma}_{R}^{2}$, where ${\sigma}_{R}$ is the median of the standard deviations ${\sigma}_{ij,R}$ of the contrast histograms ${H}_{ij}(R(x))$. The neighborhood $C(x)$ is defined as x itself and the four 26connected neighboring voxels with highest probability $P(y\sim xR)$, thus representing the neighbors most likely to be connected to x. The model is similar to a diffusion process and can be estimated with an iterated conditional modes (ICM) approach, updating sequentially the probabilities (Bazin and Pham, 2007):
from the initial voxelwise posteriors until the ratio of changed parcellation labels decreases below 0.1, typically within 50–80 iterations. The diffused probabilities and parcellation are shown in Figure 11B.
Topology correction
Request a detailed protocolThe final step of the parcellation algorithm takes a global view of the individual structures, growing from the highest posterior values inside toward the boundaries. This region growing approach makes the implicit assumption that posterior maps should be monotonically decreasing from inside to outside, which is not necessarily the case. Therefore, we perform first a topology correction step on the individual structure posteriors $P(x\in iR,B,S,C)={\mathrm{max}}_{ij}P(x\in {B}_{ij}R,S,C)$ with a fast marching algorithm (Bazin and Pham, 2007). While the corrected posterior is very similar to the original one (see Figure 11C), it ensures that all regions obtained by growing to a threshold have spherical object topology.
Anatomical region growing
Request a detailed protocolLast, we turn the posteriors into optimized parcellations, by growing them concurrently (to avoid overlaps) until the target volume for each structure is reached. Given the volume ${V}_{i}(R,B,S,C)$ of the parcellation of the diffused and topologycorrected posteriors, we define the following target volume:
taking a weighted average of the volume estimated from the data and the prior volume. This approach ensures that even in extreme cases where some structures have low posteriors, they are still able to grow to a plausible size. The region growing algorithm is driven from the most likely voxels, defined as $P(x\in iR,B,S,C){\mathrm{max}}_{j\ne i}P(x\in jR,B,S,C)$, and further modulated to follow isocontours of the skeleton prior:
Directionality of internal structures is a useful tool for understanding mechanical function in bones (Maquer et al., 2015). Here, we adapt this concept by using the skeleton isocontours as a representation of internal directionality, maintaining the intrinsic shape of structures. Thus, voxels with highest probability compared to the other structures and with similar distance to the internal skeleton are preferentially selected. The final parcellation is given in Figure 11C.
Validation metrics
Request a detailed protocolTo validate the method against manual expert delineations, we compared the MASSP results and the expert delineations with the following three measures:
The Dice overlap coefficient (Dice, 1945) $D(A,B)=\frac{2A\cap B}{A+B}$, which measures the strict overlap between voxels in both delineation;
The dilated overlap coefficient $dD(A,B)=\frac{A\cup d(B)+B\cup d(A)}{A+B}$, where $d(.)$ is a dilation of the delineation by one voxel, which measures the overlap between delineations allowing for one voxel of uncertainty;
The average surface distance $asd(A,B)$, measuring the average distance between voxels on the surface boundary of the first delineation to the other one and reciprocally, which measures the distance between both delineations.
We computed all three measures for the manual delineations from the two independent raters, as well as the ratio of overlaps (automated over manual) and distances (manual over automated) to compare both performances, as detailed in the Results section.
Comparisons with other automated methods
Request a detailed protocolTo assess the performance of MASSP compared to existing parcellation tools, we ran Freesurfer (Fischl et al., 2002), FSL FIRST (Patenaude et al., 2011) and a multiatlas registration approach (coregistering 9 of the 10 manually delineated subjects on the remaining one with ANTs [Avants et al., 2008] and labeling each structure by majority voting, similarly to the MAGeT Brain approach of Chakravarty et al., 2013). Freesurfer and FIRST were run on the skullstripped R1 map, while the multiatlas approach used all three R1, R2*, and QSM contrasts. All methods were compared in terms of Dice overlap, dilated overlap and average surface distance. We also assessed the presence of a systematic volume bias, defined as the average of the signed difference of the estimated structure volume to the manually delineated volume, normalized by the manually delineated volume.
Application to new MRI contrasts
Request a detailed protocolBefore applying MASSP to unseen contrasts, we need to convert its intensity prior histograms ${H}_{ij,R}$ to the new intensities. In order to perform this mapping, we first created a groupwise median of the HCP subjects, by coregistering every subject to the MASSP template using ANTs with nonlinear registration and both T1w, T2w contrasts matched to the template’s R1 and R2* maps. The histogram bins are then updated as follows:
adding the joint probability of the quantitative contrasts weighted by their importance for each interface to define the new intensity histograms. This model is essentially projecting the joint likelihood of the MASSP contrasts onto the new contrasts, assuming that the coregistration between the two is accurate enough. With these new histograms, we compared the testretest reliability and overall agreement of MASSP with Freesurfer parcellations included in the HCP preprocessed data set.
Measurement of structure thickness
Request a detailed protocolFinally, when comparing derived measures obtained over the lifespan with MASSP compared to manual delineations, we explored the utility of a shape thickness metric, based on the medial representation. Given the signed distance function ${\phi}_{i}$ of the structure boundary and s_{i} of the structure skeleton, the thickness is given by:
Like in cortical morphometry, thickness is a local measure, defined everywhere inside the structure, and expected to provide additional information about anatomical variations. Indeed, a similar measure of shape thickness has recently been able to highlight subtle anatomical changes in depression (Ho et al., 2020).
Software implementation
Request a detailed protocolThe proposed method, Multicontrast Anatomical Subcortical Structure Parcellation (MASSP), has been implemented as part of the Nighres toolbox (Huntenburg et al., 2018), using Python and Java for optimized processing. The software is available in open source from (release1.3.0) and . A complete parcellation pipeline is included with the Nighres examples. Computations take under 30 min per subject on a modern workstation.
Data availability
The tool presented in this article is available in open source on Github (https://github.com/nighres/nighres). The atlases necessary to run the algorithm have been deposited on the University of Amsterdam FigShare (https://doi.org/10.21942/uva.12074175.v1 and https://doi.org/10.21942/uva.12301106.v2). A single sample subject data set has been deposited on the University of Amsterdam FigShare (https://doi.org/10.21942/uva.12280316.v2). All the measurements used to generate the figures included in the article have been deposited on the University of Amsterdam FigShare (https://doi.org/10.21942/uva.12452444.v1).

FigShare / University of Amsterdam / Amsterdam University of Applied SciencesThe Amsterdam Ultrahigh field adult lifespan database (AHEAD): A freely available multimodal 7 Tesla submillimeter magnetic resonance imaging database.https://doi.org/10.21942/uva.10007840.v1

WUMinn HCP Retest DataID hcpyoungadult. The WUMinn Human Connectome Project.
References

Iron deposits in the chronically inflamed central nervous system and contributes to neurodegenerationCellular and Molecular Life Sciences 71:1607–1622.https://doi.org/10.1007/s0001801315098

BookDigital Homeomorphisms in Deformable RegistrationIn: Karssemeijer N, Lelieveldt B, editors. Information Processing in Medical Imaging. Berlin, Heidelberg: Springer Berlin Heidelberg. pp. 211–222.https://doi.org/10.1007/9783540732730_18

ConferenceVessel segmentation from quantitative susceptibility maps for local oxygenation venography2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI). pp. 1135–1138.https://doi.org/10.1109/ISBI.2016.7493466

Denoising HighField MultiDimensional MRI with local complex PCAFrontiers in Neuroscience 13:66.https://doi.org/10.3389/fnins.2019.01066

Topology correction of segmented medical images using a fast marching algorithmComputer Methods and Programs in Biomedicine 88:182–190.https://doi.org/10.1016/j.cmpb.2007.08.006

Homeomorphic brain image segmentation with topological and statistical atlasesMedical Image Analysis 12:616–625.https://doi.org/10.1016/j.media.2008.06.008

Biological shape and visual science. IJournal of Theoretical Biology 38:205–287.https://doi.org/10.1016/00225193(73)901756

MP2RAGEME: T 1 , T 2 * , and QSM mapping in one sequence at 7 teslaHuman Brain Mapping 40:1786–1798.https://doi.org/10.1002/hbm.24490

A unified framework for crossmodality multiatlas segmentation of brain MRIMedical Image Analysis 17:1181–1191.https://doi.org/10.1016/j.media.2013.08.001

Critical ages in the life course of the adult brain: nonlinear subcortical agingNeurobiology of Aging 34:2239–2247.https://doi.org/10.1016/j.neurobiolaging.2013.04.006

Towards a mechanistic understanding of the human subcortexNature Reviews Neuroscience 18:57–65.https://doi.org/10.1038/nrn.2016.163

Ultra high field MRIGuided deep brain stimulationTrends in Biotechnology 35:904–907.https://doi.org/10.1016/j.tibtech.2017.06.010

Structural brain magnetic resonance imaging of limbic and thalamic volumes in pediatric bipolar disorderAmerican Journal of Psychiatry 162:1256–1265.https://doi.org/10.1176/appi.ajp.162.7.1256

Ultrahigh 7T MRI of structural agerelated changes of the subthalamic nucleusJournal of Neuroscience 33:4896–4900.https://doi.org/10.1523/JNEUROSCI.324112.2013

Effects of aging on T₁, T₂*, and QSM MRI values in the subcortexBrain Structure & Function 222:2487–2505.https://doi.org/10.1007/s0042901613524

Bone volume fraction and fabric anisotropy are better determinants of trabecular bone stiffness than other morphological variablesJournal of Bone and Mineral Research 30:1000–1008.https://doi.org/10.1002/jbmr.2437

Atlas generation for subcortical and ventricular structures with its applications in shape analysisIEEE Transactions on Image Processing 19:1539–1547.https://doi.org/10.1109/TIP.2010.2042099

Fast robust automated brain extractionHuman Brain Mapping 17:143–155.https://doi.org/10.1002/hbm.10062

BookHandbook of Basal Ganglia Structure and FunctionIn: Steiner H, Tseng K, editors. Handbook of Behavioral Neuroscience, 24. Elsevier. pp. 1–1036.

Myelin and modeling: bootstrapping cortical microcircuitsFrontiers in Neural Circuits 13:34.https://doi.org/10.3389/fncir.2019.00034

ConferenceGeneration and evaluation of an ultrahighfield atlas with applications in DBS planningProceedings of SPIE Medical Imaging.

Multicontrast unbiased MRI atlas of a parkinson's disease populationInternational Journal of Computer Assisted Radiology and Surgery 10:329–341.https://doi.org/10.1007/s115480141068y

Lifespan maturation and degeneration of human brain white matterNature Communications 5:4932.https://doi.org/10.1038/ncomms5932
Article and author information
Author details
Funding
Nederlandse Organisatie voor Wetenschappelijk Onderzoek (VICI)
 Birte U Forstmann
Nederlandse Organisatie voor Wetenschappelijk Onderzoek (STW)
 Anneke Alkemade
 Birte U Forstmann
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Josephine Groot, Nikita Berendonk, Nicky Lute for their help collecting the AHEAD database, and Wietske van der Zwaag and Matthan Caan for their help in setting up the MP2RAGEME sequence. We also thank Steven Miletić and Dagmar Timmann for stimulating discussions around this topic, and all undergraduate students who contributed to the manual delineations. This work was supported by a NWO Vici grant (BF) and a NWO STW grant (AA, BF). HCP data were provided by the Human Connectome Project, WUMinn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.
Ethics
Human subjects: Informed consent and consent to publish, including consent to publish anonymized imaging data, was obtained for all subjects. Ethical approval was obtained with the University of Amsterdam Faculty of Social and Behavioral Sciences LAB Ethics Review Board, with ERB number 2016DP6897.
Version history
 Received: May 28, 2020
 Accepted: December 15, 2020
 Accepted Manuscript published: December 16, 2020 (version 1)
 Version of Record published: December 29, 2020 (version 2)
Copyright
© 2020, Bazin et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,433
 views

 227
 downloads

 24
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
 Neuroscience
Normal aging leads to myelin alterations in the rhesus monkey dorsolateral prefrontal cortex (dlPFC), which are positively correlated with degree of cognitive impairment. It is hypothesized that remyelination with shorter and thinner myelin sheaths partially compensates for myelin degradation, but computational modeling has not yet explored these two phenomena together systematically. Here, we used a twopronged modeling approach to determine how agerelated myelin changes affect a core cognitive function: spatial working memory. First, we built a multicompartment pyramidal neuron model fit to monkey dlPFC empirical data, with an axon including myelinated segments having paranodes, juxtaparanodes, internodes, and tight junctions. This model was used to quantify conduction velocity (CV) changes and action potential (AP) failures after demyelination and subsequent remyelination. Next, we incorporated the single neuron results into a spiking neural network model of working memory. While complete remyelination nearly recovered axonal transmission and network function to unperturbed levels, our models predict that biologically plausible levels of myelin dystrophy, if uncompensated by other factors, can account for substantial working memory impairment with aging. The present computational study unites empirical data from ultrastructure up to behavior during normal aging, and has broader implications for many demyelinating conditions, such as multiple sclerosis or schizophrenia.

 Neuroscience
How is new information organized in memory? According to latent state theories, this is determined by the level of surprise, or prediction error, generated by the new information: a small prediction error leads to the updating of existing memory, large prediction error leads to encoding of a new memory. We tested this idea using a protocol in which rats were first conditioned to fear a stimulus paired with shock. The stimulus was then gradually extinguished by progressively reducing the shock intensity until the stimulus was presented alone. Consistent with latent state theories, this gradual extinction protocol (small prediction errors) was better than standard extinction (large prediction errors) in producing longterm suppression of fear responses, and the benefit of gradual extinction was due to updating of the conditioning memory with information about extinction. Thus, prediction error determines how new information is organized in memory, and latent state theories adequately describe the ways in which this occurs.