Atlas-guided segmentation of vervet monkey brain MRI.

The vervet monkey is an important nonhuman primate model that allows the study of isolated environmental factors in a controlled environment. Analysis of monkey MRI often suffers from lower quality images compared with human MRI because clinical equipment is typically used to image the smaller monkey brain and higher spatial resolution is required. This, together with the anatomical differences of the monkey brains, complicates the use of neuroimage analysis pipelines tuned for human MRI analysis. In this paper we developed an open source image analysis framework based on the tools available within the 3D Slicer software to support a biological study that investigates the effect of chronic ethanol exposure on brain morphometry in a longitudinally followed population of male vervets. We first developed a computerized atlas of vervet monkey brain MRI, which was used to encode the typical appearance of the individual brain structures in MRI and their spatial distribution. The atlas was then used as a spatial prior during automatic segmentation to process two longitudinal scans per subject. Our evaluation confirms the consistency and reliability of the automatic segmentation. The comparison of atlas construction strategies reveals that the use of a population-specific atlas leads to improved accuracy of the segmentation for subcortical brain structures. The contribution of this work is twofold. First, we describe an image processing workflow specifically tuned towards the analysis of vervet MRI that consists solely of the open source software tools. Second, we develop a digital atlas of vervet monkey brain MRIs to enable similar studies that rely on the vervet model.

INTRODUCTION

Nonhuman primate (NHP) models are essential in understanding brain function and disease. The macaque and African green (vervet) monkeys are common NHP models in translational research. Monkeys are close to humans physiologically and behaviorally, justifying their use in studies that cannot be done on humans. Additionally, control of their experimental history and environment avoids many of the variables that confound clinical studies. The vervet monkey is one of the most important NHP models, which is becoming increasingly used [1-5]. Compared to rhesus macaques, vervets are abundant, less expensive, and are easy to handle [6]. Nonhuman primates [7], including the vervet monkey [8-10], have been increasingly used to provide a unique animal model for studying the effects of alcohol abuse. The vervet monkey (Chlorocebus aethiops) is anecdotally thought to consume ethanol in nature and as such, provides a relevant model of the neurobiological consequences of long-term ethanol exposure.Our work is motivated by a biological study focused on the effects of long-term heavy ethanol (EtOH) consumption on the vervet brain structure using in vivo Magnetic Resonance Imaging (MRI) [11]. In order to quantify these effects, we are interested in measuring the morphological changes of the individual brain structures over time.

In the context of the biological study, we have been following a population of 10 vervets, which include normal controls that consume a maltose-dextrin solution that is isocaloric to ethanol and those subjects that are induced to consume ethanol, as previously described by Vivian et al. [12]. The subjects are imaged prior to ethanol or maltose-dextrin induction (baseline imaging), and at a post-induction interval (follow-up imaging). The ultimate objective of our image analysis is to provide accurate volumetric measurements of certain anatomical structures thought to be impacted by chronic alcohol ingestion. Morphological analysis has been used extensively for human brain MRI [13,14]. However, the existing tools are fine-tuned to human brain MRIs, which are quite different from monkey brain MRIs. Compared to human brain MRIs, scans of the monkey brain are acquired with higher spatial resolution to capture anatomical detail of the smaller brain. This often results in a poorer signal-to-noise ratio (SNR), as compared to human MRIs. The imaging coils used for NHP MRIs are typically not developed specifically for NHP imaging (as an example, human knee coils are sometime used to image the NHP brain [15]) or are not as advanced as those used for humans. This contributes to considerable signal inhomogeneity and can reduce SNR (see Fig. ). Most importantly, there are notable anatomical differences between the monkey and human brain. Specifically, the olfactory cortex of the NHP is more prominent at rostral levels than in the human brain, NHP have smaller size of the cerebral cortex that has different gyration patterns, and the relative size of the subcortical structures is different between the human and NHP brains [16]. A central issue in adapting the existing tools to the task of processing NHP MRI is the generation of computerized atlases of the brain anatomy [17], which are used to guide automated segmentation tools. These atlases encode different characteristics of the brain anatomy for the studied population, such as the MR signal intensity, shape and the relative location of the individual brain structures with respect to each other [17,18]. This information is especially critical for automatic segmentation of the anatomical structures with similar MR intensity signal profiles. While atlases for some NHP species have been developed previously [15,19-24], no such atlas was available for the vervet at the time of our study. In this paper we develop such an atlas, and integrate it into an automatic atlas-guided segmentation tool, EM Segmenter [25] available within 3D Slicer [26,27].

Constructing computerized MRI atlases often requires segmentation of a representative set of scans (atlas population), and aligning them to some reference. The spatial distribution of tissue classes is then computed by calculating the frequency of occurrence for each tissue class in the spatially aligned segmentations. Manually segmenting the MRIs of our vervet population was particularly difficult due to the low signal to noise ratio and strong bias artifacts. These artifacts cause manual segmentation of individual scans to be time-consuming, and to have a high inter- and intra-rater variability. We address this issue by following an alternative strategy originally suggested by Styner et al. [19,28]. We constructed an average segmentation template image by spatially aligning the baseline scans of the atlas population to a common reference and averaging them. The resulting average image has reduced noise, which simplifies and increases the reproducibility of manual segmentation, making it more practical. The transformations that align the individual scans with the reference image are used to back-propagate the manual segmentation to each of the subjects in the atlas population. The probabilistic atlas is constructed by computing frequencies of occurrence for the individual tissues in these back-propagated segmentations. We use two different strategies to align the individual images. The first method uses a fixed reference subject, while the second approach uses an unbiased strategy [29-31], where the template is iteratively updated so that each subject in the atlas population contributes equally to the atlas.

In order to evaluate our segmentation workflow, we apply it to the baseline images, and to the scans acquired during the first follow-up imaging session. Our focus is on the development and evaluation of the segmentation technology, therefore in this paper we do not attempt to quantify the structural changes caused by the exposure of subjects to alcohol. However, since the follow-up images were not used for atlas construction, they enable us to evaluate the robustness of the workflow while processing images that have different properties, and that reflect longitudinal changes in the subject anatomy.

In summary, we present a pipeline consisting of open source tools available within 3D Slicer, which were carefully tuned towards segmenting vervet brain MRIs. We evaluate this pipeline on two sets of MRIs for our vervet population, and study the impact of different atlas construction methods on the segmentation result. Finally, we have published online the digital atlas of the vervet monkey brain [32] and provide in this article a detailed description of the pipeline itself. This hopefully will reduce the technical barriers in automatically segmenting vervet brain MRIs.

MATERIALS AND METHODS

Our study uses two sets of MRIs, which were acquired before and after induction of the animals to self-administration of EtOH or an isocaloric maltose-dextrin solution, as described in Section 2.1. Section 2.2 presents our approach to automatic segmentation of the vervet brain. First, we preprocess the raw baseline and follow-up MR images to improve the reliability of the subsequent analysis. The preprocessed baseline scans are used for the construction of population-specific digital atlas capturing the typical appearance of the normal vervet brain in MRIs. The details of integrating this atlas into an automatic segmentation tool within 3D Slicer are discussed in Section 2.2. Finally, we describe our evaluation methodology and the statistical methods in Section 2.3.

Image Data

Ten male vervet monkeys (Chlorocebus aethiops) were obtained from the Wake Forest University Primate Center (Winston-Salem, NC) and singly housed in cages equipped with operant panels through which all food and fluids were provided. Baseline imaging was performed when all animals were approximately 4 years old. All animals were trained to operate drinking panels in an identical fashion in an ongoing ethanol self-administration study funded by National Institute on Alcohol Abuse and Alcoholism (NIAAA) to Dr. David Friedman. The animals were trained to self-administer either EtOH, or an isocaloric maltose-dextrin (control) solution through the panel. All procedures in this study were conducted in compliance with State and Federal laws, standards of the US Department of Health and Human Services, and guidelines established by the Wake Forest University Institutional Animal Care and Use Committee as well as the National Institute of Health Guide for the Care and Use of Laboratory Animals (NIH Publications No. 80-23)

The imaging data was collected under conditions of general anesthesia. Briefly, all animals were initially sedated with ketamine and transported to the MRI Center where they were placed on the scanner bed. They inhaled isoflurane (3%) through a nose-cone at which time they were then intubated and the concentration of isoflurane was reduced to 1.5% for the remainder of the scans. The animals were artificially ventilated to maintain the same physiological range of vital signs across all animals. Vital signs including oxygenation saturation, heart rate, respiration rate, expired CO2 and isoflurane concentration were monitored every 10-15 minutes using a pulseoximeter and an anesthesia monitor for the duration of the scans. At the end of the scan session, the isoflurane flow and ventilator were shut off and the animals were allowed to breathe a mixture of oxygen and room air until they began to recover from the anesthesia. Once they exhibited various signs of recovery such as tongue and eye movements and gag reflex, they were extubated and provided pipeline oxygen until they recovered enough to sit in an upright position. They were then placed back into an animal carrier and returned to their home cages where they were further monitored until they were fully alert.

Axial T1-weighted MR images were acquired on a 3T GE scanner with a circularly-polarized, single channel dedicated RF coil with an internal diameter of 18.4 cm (Litzcage, Doty Scientific, Columbia, SC), using a 3D SPGR sequence (TI 600ms, TE 3.276ms, TR 15.28ms; flip angle 15o; matrix 256x256; FOV 12cm; in-plane resolution 0.47mm; slice thickness 0.5mm). The raw images were reconstructed into 3D volumes before processing. Two sets of scans were collected. The baseline scan was acquired at the age of 4 years. The follow-up acquisition was done approximately 8 months later. During the first 4 months between the baseline and follow-up acquisition the monkeys were trained to operate the drinking panel, followed by the 4 month induction period [33]. Only half of the studied vervet population was exposed to ethanol and underwent the EtOH induction procedure. The other half was the control group that operated the drinking panel in an identical fashion as the EtOH group but consumed an isocaloric maltose-dextrin solution.

Image Segmentation Pipeline

Our image processing is based on the software components implemented in 3D Slicer [27] and ANTS [34]. 3D Slicer is a free open source software for medical image computing and visualization. The image processing needs within this project included visualization and annotation of the images, ability to process DICOM data (needed for the conversion of the raw data obtained from the scanner), interface for manual image segmentation (preparation of the manually segmented template image), linear and deformable image registration capabilities that could be customized to our data, automatic segmentation of the brain MRIs, and the ability to develop new image processing components that can be integrated within the framework. 3D Slicer provides all of the needed functionality within the unified application interface.

Our segmentation pipeline consists of the three main components: preprocessing, atlas formation, and segmentation. During preprocessing we perform spatial alignment of the images, skull stripping, and correction for intensity inhomogeneity. Next we prepare a population-specific atlas of vervet brain MRI constructed from the baseline scans. As part of atlas construction, we compute the spatial priors for selected brain regions, which will guide the segmentation process. We consider two approaches to atlas construction. In the first case, the biased template is prepared by spatially registering all of the images in the baseline population to one of the subjects. In the second case, the template is derived by iterative registration of the baseline scans to their pixelwise average, resulting in an unbiased atlas template of our population. We then register this template to each of the baseline and follow-up scans, and apply the resulting transformation to align the probabilistic atlases. The transformed atlases serve as a spatial prior in the final step to guide the EM Segmenter, which we specifically parametrized to NHP MRIs.

Image Preprocessing

Our preprocessing of the raw MR images included skull-stripping, correction of the bias field inhomogeneity artifact and correction for the head pose, as summarized in Fig. (). Despite the continuous advancements in the development of skull-stripping methodology [35-37], there are few robust tools that are widely available. The available skull-strippers, e.g., those available in FSL and SPM [13,14], are fine-tuned for processing human brain anatomy, and the results are highly dependent on the anatomical features of a subject [28]. Several existing studies applied BET, the skull-stripping tool distributed as part of FSL, to brain MRIs of NHP [15,22,24]. However, due to the differences in the head anatomy in human and NHP, this requires careful re-parameterization of BET. In a recent study, BET was applied to the challenging task of skull-stripping in the MRI of baboon fetus [24]. The authors note that even after careful parameterization, “manual editing is often necessary to insure anatomically correct brain extraction”. Our preliminary evaluation of applying BET to vervet brain MRI resulted in a similar conclusion, which motivated us to develop a customized template-based approach to skull-stripping using the tools available within 3D Slicer, as shown in Fig. ().

A representative vervet brain MRI (ICC registration template) was manually skull-stripped (Step 1 of the pre-processing workflow, see Fig. ), and then non-rigidly registered to each of the remaining MRIs. The obtained transformation was applied to the ICC mask in the template, providing a good initial approximation of the ICC in the individual subject, which was then manually edited if necessary. A critical component of this approach is inter-subject non-rigid registration. While the appearance of the vervet brain is not significantly different across subjects, this is not the case for the extra-cranial structures and tissues. The thickness of the temporalis muscle equals or exceeds the width of the brain at more rostral locations, and is approximately the same as the brain width at more caudal locations. Appearance and shape of the soft extra-cranial tissues (specifically, the temporalis muscle, neck muscles, ears and the jaw) are usually not consistent between the image acquisitions. Perhaps more important are the differences in the location of the field of view relative to the subject head, which lead to variable portions of the neck being present in the image. The vervet brain occupies a relatively small portion of the image volume (~5% of the image voxels), therefore in the naïve registration between the two vervet MRIs (using all of the image voxels), those voxels that correspond to the brain contribute little weight to the overall value of the similarity metric.

In order to achieve robust registration between the ICC template and the reference MRI (the subject being skull-stripped), the following customizations of registration proved useful. We automatically segment the head region in the reference image using BRAINSROIAuto tool of 3D Slicer (Step 2 in Fig. ). The N4ITK algorithm [38] is applied to both the ICC template and the reference MRI to compensate for the acquisition-specific spatial differences in the signal intensity in the brain region (Step 3). The head region mask is used by N4ITK to improve the robustness of the inhomogeneity correction. The bias-corrected images are next automatically registered using the BRAINSFit module [39] of 3D Slicer (Step 4). The registration strategy implemented in BRAINSFit searches for a transformation that minimizes the negated Mutual Information image similarity measure [40] using a gradient descent optimizer. The registration process is parameterized by the ICC mask in the template image, and the head mask in the reference image: the calculation of the similarity metric is restricted to the overlap of the regions defined by these two masks. The complexity of the deformation model is gradually increased from rigid, to affine, to B-spline transformations. In some cases initial rigid aligment of the template to the reference image was necessary due to the large discrepancy in the relative head location or orientation. This was accomplished using the interactive registration tools available in 3D Slicer. Following the registration of the template to the reference image, the ICC template mask was resampled using the final transformation. The result was examined visually, and minor manual corrections were performed using the Editor module of 3D Slicer to insure accurate skull-stripping (Step 4).In the last pre-processing step (Step 5) we positioned all images into consistent orientation. The ICC template image was manually oriented so that the AC-PC line was orthogonal to the coronal plane, and the midline plane was aligned with the sagittal plane. Three landmarks corresponding to the AC, PC and a point on the midline were manually identified in each of the images. A closed-form solution [41] for rigid transformation that aligns these three landmarks in each of the subjects with landmarks in the ICC template was automatically calculated using 3D Slicer. Finally, we rescaled the intensities in each of the skull-stripped images to have the same range.

As the result of this preprocessing, we obtained skull-stripped, bias-corrected and consistently oriented scans across baseline (pre-induction) and follow-up (post-induction) MRIs of the vervet population.

Atlas Construction

We use the pre-processed scans to capture the typical anatomy of the vervet brain in a population-specific digital atlas. This atlas provides the prior information about the location of individual brain structures to the EM segmentation algorithm of Pohl et al. [25]. The atlas consists of two components. A probabilistic atlas defines the probability of occurrence of a given tissue for each voxel of the image. It is defined in the space of the average anatomy of the studied population, represented by an average intensity atlas (often referred to as atlas template) that stores the typical signal produced by a given tissue during MRI acquisition. In our workflow, the atlas is constructed from the single manually segmented template. We use the non-rigid transformations between scans in the atlas population and this template to estimate normal variability of the brain, as we discus below. In addition to the two population-specific atlases, we use the rhesus atlas developed by Styner et al. [19] for comparison. Since the vervet is similar anatomically to the rhesus monkey, it is a reasonable candidate atlas for segmenting vervet brain MRIs. We show representative slices from the three atlases used in this study in Fig. (), and axial slices at different locations for the unbiased population-specific atlas in Fig. (). Our atlas population is restricted to scans obtained before the induction procedure. This ensures that we do not encode atypical appearance of the vervet brain that may be caused by exposure to alcohol.

Usually, atlas construction requires several manually segmented images to derive the statistics of spatial distribution of tissues [25]. Manual segmentation is a tedious and error-prone task. Based on our experience, manual segmentation of a single structure in one image requires on average 30 minutes of operator time, and is complicated by inconsistent image quiality and noise. Significant differences were observed in the contrast-to-noise ratio between different subjects and between the different imaging sessions due to the continuous improvements of the imaging protocols throughout the study, as shown in Table . An alternative approach, which does not require all of the atlas population scans to be segmented, was suggested by Styner et al. [19]. Using this approach, segmentation is prepared manually or semi-automatically for a single scan. The result is back-propagated onto each subject in the atlas population using non-rigid registration. These back-propagated segmentations are used to construct the probabilistic atlas in the same fashion as in the conventional atlas construction approaches. Next we describe this atlas construction procedure in more detail.

The baseline scan that was least affected by the imaging artifacts and had high signal to noise ratio was selected as segmentation template. This was done to facilitate manual segmentation of this scan. The remaining scans of the baseline population were registered to this reference volume using BRAINSFit, as described in the previous section. The transformation recovered by the registration for each subject consisted of the two components: an affine and a B-spline transformation. The affine component of the transformation encodes the global differences between the template and each of the subjects (e.g., due to the difference in pose and total brain volume), while the non-rigid component represents the subject-specific differences in individual brain structures (e.g., differences in the cortex gyration). The average template was computed by voxel-wise averaging of the spatially aligned T1w MRI images. The template segmentation was initialized using k-means clustering available in Insight Toolkit [42] into three classes corresponding to cerebrospinal fluid (CSF), grey matter (GM) and white matter (WM). This initial segmentation was examined by a trained operator and manually corrected as necessary using the image editing capabilities of 3D Slicer. The subcortical structures of interest (putamen (PUT), caudate (CAU) and hippocampus (HPC)) were traced manually in one hemisphere. Segmentations of subcortical structures in the other hemisphere were obtained by mirroring across the midsaggital plane and adjusting the result manually.

Next, we computed the inverse of the B-spline non-rigid component of the transformation for each of the subjects, which were previously registered to a single subject in the population, using the iterative inverse deformation field estimation method implemented in Insight Toolkit [42]. The resulting transformation was used to back-propagate segmentation of the average template to each of the subjects in the atlas population. Voxel-wise probability maps for each of these tissues were obtained by counting their frequencies of occurrence in the back-propagated segmentations. This approach aligns all of the scans in the atlas population to the single subject. The resulting atlas is biased towards the anatomy of the selected reference scan: the average atlas template will be more similar to the reference subject compared to other subjects in the population. To address this issue we also generated an atlas using the unbiased template generator developed by Avants and Gee [31] (available in the open source ANTS library [34]). This approach iterates between the computation of the linear average of the scans in the baseline population, and diffeomorphic registration of each of the scans to this average. We applied this method to the set of baseline scans, and obtained an unbiased average template, a deformation field that maps each of the scans in the atlas population to this average, and an inverse of this deformation field. The probabilistic atlas was constructed using back-propagation as described earlier. First, we warped the manual segmentation prepared for the biased template to the unbiased template. Next we applied the inverse transformations produced by the unbiased atlas estimation to back-propagate the segmentation from the unbiased template to each of the scans in the atlas population. The probabilistic atlas was computed by calculating the frequencies of occurrence for each of the segmentation labels.

As the final step of atlas processing, we non-rigidly registered the average atlas template to each of the baseline and follow-up scans. The transformation recovered by the registration was applied to warp the corresponding probabilistic atlases for each of the tissue classes to the anatomy of the subject. These aligned probabilistic atlases were used in the next step of the processing pipeline to guide automatic segmentation.

Configuration of the Segmentation Algorithm

The final step of our processing pipeline is automatic segmentation of the images. We use the EM segmentation approach by Pohl et al. [25] as it enables relatively easy adaptation for our task. This approach is parameterized using a rooted tree data structure. Using graph theory terminology, tree is defined as a connected graph without cycles. The graph edges in the rooted tree used by EM Segmenter are directed away from the root node that describes the region being segmented (e.g., the brain). Tree leaves (the nodes that do not have outgoing edges) are located at the bottom of the hiearchy. Intermediate tree nodes in between the root and leaf nodes correspond to the logical groups of the anatomical structures (e.g., left subcortical gray matter), while the leaves represent the individual segmented structures (e.g., hippocampus). In addition to encoding the parent-child relationships between the segmented structures, the tree stores priors for each individual tissue class. For example, segmentation of the hippocampus proceeds by separating the brain volume into gray matter, white matter and cerebrospinal fluid classes, segmenting the subcortical structures from the total gray matter volume, and finally segmenting the hippocampus from the subcortical gray matter class. Such organization of the segmentation process improves the robustness and simplifies parameterization. The configuration of the segmentation hierarchy tree was customized for the purposes of segmenting the population of vervet MRI scans using the hierarchy shown in Fig. (). The tree was populated with the subject-specific information: a spatial prior, which was prepared in the previous step, and an intensity model. The approximate intensity model (mean and variance of the intensity distribution for the given tissue class) was estimated automatically by sampling the intensities for each structure in the areas where it is likely to be located (95% probability, as defined by the spatial prior).

In addition to the spatial prior and intensity model, segmentation of each tissue class was guided by the confidence parameter λ in the range zero to one. At each leaf of the tree the confidence specifies whether the classification decision should rely more on the spatial prior (λ=1), or on the intensity model (λ=0). The idea behind the specific organization of the hierarchical tree (see Fig. ) was to use similar values of the confidence parameter for each structure at a given hierarchy level, thus simplifying parameterization of the algorithm. On levels 2, 3, 5 and 6 we set λ value close to 1 (intensities are very similar for brainstem and white matter, thus we need to rely on the spatial prior). Segmentation of white/gray matter and CSF boundaries relied more on the intensity distribution model, thus we set λ≈0 at this level. The specific values of parameters were derived by experimentation on a single image from the baseline population, and then used without any additional tuning to segment the rest of the images. We implemented the described configuration of the segmentation algorithm in a customized NHP segmentation task of the EM Segmenter module of 3D Slicer [43].

Evaluation Methodology

In our evaluation we studied the agreement between the automatic segmentations and the human raters, and the effect of the atlas choice on this agreement. The purpose of the evaluation of the workflow on the follow-up MRIs was to confirm that it can be applied to the scans not used during the atlas construction.

The quantitative comparison used three rectangular regions (further referred as R1, R2 and R3) defined in the coronal plane, manually segmented in each image by two trained research assistants using 3D Slicer, see Fig. (). The regions include sections of each of the subcortical structures that were of interest for the biological study (caudate, hippocampus and putamen). R1 was in the frontal part of the brain just anterior to the lateral ventricles. The bottom left corner of R2 was at the anterior commissure, capturing a cross-section of the putamen and caudate nucleus. R3 was located at the coronal slice through posterior commissure, intersecting hippocampus and temporal lobe. In each of the regions the boundaries of GM, WM, putamen, caudate and hippocampus were manually traced.

Our evaluation of agreement between the results of automatic and manual segmentation used the Dice Similarity Coefficient (DSC) [44] and Multidimensional scaling (MDS) technique. DSC provides a numeric measure of agreement in the range between 0 and 1, with 1 corresponding to the perfect agreement. DSC can be used to quantify the agreement between two labels. Multidimensional scaling (MDS) [45] can be used in combination with the pairwise DSC assessment to explore the relative performance of more than two segmentation methods. MDS is a mathematical technique, which finds the optimal projection of a collection of points in a high-dimensional space into two dimensions in such a way that the distances between corresponding projected points of the graph are most similar to the distances in the original higher-dimensional space [46]. Applied to our evaluation, the set of points used in MDS analysis corresponded to the five segmentation results of the validation regions R1-R3 (three obtained using automatic segmentation with different atlases, and two completed by human raters). All three segmented regions corresponding to the individual subject were merged into one volume for the calculation of DSC, so that the DSC measure for white matter, gray matter and CSF was computed over the labels of these structures from all three regions. The distance between the compared segmentation results was calculated for each individual structure as the dissimilarity (1-DSC) averaged across the 20 scans of the subjects (pooled baseline and follow-up) for each of the labels. Multidimensional scaling was performed using standard functionality of Matlab 7.10 (Mathworks, Natick, MA, USA). The accuracy of the MDS projections was examined using the residual plots that show the correspondence between the original and projected distances between the points. In the cases when the projection is good, the projected and original distances between the points should be similar.

The MDS projections provide a summary view of the agreement between multiple pairs of measurements. For the purposes of quantitative evaluation, we compared the DSC measure of overlap between the results of the manual tracing and the segmentations obtained automatically for gray matter, white matter, putamen, caudate and hippocampus. Given the DSC measures computed between the segmentations produced by one of the raters and each of the automatic results obtained using different atlases, we performed group comparison to quantify the differences in their agreement with the manual segmentation. We also evaluated whether there is significant difference between the segmentation accuracy for the baseline scans used to construct the atlas, and the follow-up scans.

Two-sided paired t-tests were used for evaluating the significance of the difference between two groups. One-way ANOVA, followed by the Tukey's range test, was applied to compare more than two groups. In the cases where statistical tests were applied to DSC values, the data was logit transformed [47], since the distribution of DSC does not satisfy the normality assumption. DSC measures for the baseline and follow-up scans were pooled together in this comparison. Statistical analysis was using R Environment for Statistical Computing (R Core Development Team, Vienna, Austria).

RESULTS

Three sets of results were obtained using each of the atlases considered in this study. Processing was performed on a computer equipped with a 16-core 2.6 GHz Intel Xeon CPU. Non-rigid registration of the template to the subject required approximately 11 minutes for B-spline registration, and 24 minutes for the diffeomorphic registration. Once the alignment of the spatial priors with the individual scan was recovered, the processing time required for segmenting one dataset using the multi-threaded implementation of EM Segmenter (the standard implementation available as part of 3D Slicer) was under 7 minutes. The segmentation results were visually examined to ensure overall correctness of the segmentations, see Fig. () for an example.

Representative MDS projections for the segmented tissues are provided in Fig. (). Small values of the embedding errors led us to conclude that the quality of projections was satisfactory. In most cases the results produced with the population-specific vervet atlas (using either biased or unbiased approach) were nearly equidistant from the segmentations produced by either of the manual raters. The segmentation results produced using different atlases were more consistent for gray and white matter segmentations. The segmentations of subcortical gray matter produced using the biased atlas were nearly at equal distance from each of the raters. The result obtained with the rhesus atlas for subcortical structures was more similar to the segmentations of one of the raters for the hippocampus and putamen labels. With the exception of the hippocampus, segmentations generated using the unbiased atlas were most similar to the manual segmentations, and were located at similar distance from both raters. For hippocampus, the biased atlas appears to provide results most similar to the manually produced segmentations.

Logit-transformed Dice scores for gray matter (GM), white matter (WM), putamen (PUT), caudate (CAU) and hippocampus (HPC) were compared between automatic segmentation results and manual segmentations prepared by Rater1. Overall, we observed higher variability of the agreement between the automatically generated results and the manual segmentations for the subcortical structures as can be observed in Fig. (). A two-tailed paired t-test between the DSC for the baseline and follow-up scans revealed statistically significant improvement in the overlap of the WM label segmented in the follow-up scans for all of the three atlases. This seems to indicate that the segmentation is more reliable for the follow-up scans. We believe this is due to the improved contrast-to-noise ratio in the follow-up scans (as we discussed in Section 2.2.2), which may contribute to the better quality of the automatic segmentation.

Next we used an ANOVA to assess the impact of the atlas on the agreement of the automatic segmentation result with the expert segmentation. We found that atlas choice had statistically significant effects for CAU, PUT and HPC labels. For the caudate we observed significant improvement in the DSC value due to the use of vervet-specific atlases, with the mean improvement in DSC of 0.066 (p<0.0014) and 0.083 (p<0.000015) for the biased and unbiased atlases, respectively. In the case of PUT label, we found significant improvement due to the use of the unbiased atlas instead of the biased one (mean improvement 0.063, p<0.0094). However, for the hippocampus, the only statistically significant improvement was observed due to the use of the biased instead of the rhesus atlas (mean improvement 0.057, p<0.035). The atlas choice did not have a statistically significant effect on the DSC for the WM and GM labels, in agreement with the visual assessment of the MDS analysis results. There are two contributing factors that may explain the latter observation. First, segmentation of these structures relies primarily on the signal intensity, and not on the probabilistic atlas. Second, GM and WM labels occupy larger volumes, and the DSC measure is less sensitive to the segmentation errors of the larger structures.

Our further analysis was focused on the unbiased atlas, since its use leads to improved agreement with the rater result for the caudate (compared to the rhesus atlases) and putamen labels (compared to the biased atlas), and no significant difference was found for the hippocampus. We evaluated whether the agreement between the segmentations produced automatically using the unbiased atlas and the manual tracings by each of the raters is significantly different from the inter-rater agreement. Statistical analysis across the five tissues of interest revealed that the performance of the automatic segmentation was comparable with the manually obtained results for all of the structures except the hippocampus, where the agreement between the raters is significantly higher than the agreement of the automatic segmentation results with each of the raters (one-way ANOVA followed by the Tukey's test, p<0.001).

DISCUSSION

The continuous technical improvements of MRI technology leads to increasing spatial resolution of the images, which require more time for their manual segmentation. Manual processing becomes even more infeasible as the population sizes and the number of regions of interest increases. Additionally, manual segmentation is subject to rater bias and may not be easily reproducible. The availability of image processing tools for automatic processing of MR images therefore becomes critical. In this work we narrowed the existing gap in NHP image analysis by providing an atlas of the vervet monkey brain and describing a detailed approach to vervet brain MRI segmentation using that atlas. Our evaluation showed that the performance of the automatic segmentation is comparable with the results produced by human experts. We observed particularly good quality of segmentation for the gray matter and white matter tissues, the most laborious aspect of manual tracing. We found that the automatic segmentation is less accurate for the hippocampus segmentation. The contrast between the hippocampus cross-section and surrounding tissues in R3 is stronger than for the putamen and caudate (see Fig. ), where the boundaries are diffused with the adjacent isointense gray matter. Because of this, the task of manual segmentation may be easier for this structure, which explains the better agreement between the raters for manual segmentations of the hippocampus. In the future we plan to improve automatic segmentation of this structure by incorporating shape constraints into the segmentation procedure [48].

Since the accuracy of the automatic segmentation is comparable with the results produced by the raters, the time required for the segmentation of these images can now be reduced from approximately 30 minutes of manual contouring per individual structure to under 30 minutes of unsupervised processing per scan (the time needed for the registration of the template and segmentation of the gray matter, white matter, caudate, putamen and hippocampus represented in our atlas).

We published online the average intensity atlas as well as the probabilistic atlases for the individual structures of the vervet brain [32]. The use of the population-specific atlas leads to improved agreement with the expert for all three subcortical structures, as compared to the results obtained with the atlas of a closely related rhesus species. Our processing pipeline is generic and is applicable to construction of the atlases for other species or ages, since the tools we used are publicly available. A unique aspect of our work is the quantitative comparison of the biased and unbiased atlas construction approaches. Our results show that the use of the unbiased population specific atlas results in an increased accuracy of segmentation for some of the analyzed structures, which agrees with the earlier study for human brain MRI segmentation [49].

The processing pipeline developed as the result of our study is based on open source tools, with most of the components available within 3D Slicer. We believe this will facilitate the adoption of the presented methodology in similar studies. We note that most of the existing studies in NHP MRI processing used publicly available software components. Wisco et al. [50] used 3D Slicer for some of the processing steps in their pipeline. Styner et al. [19] applied publicly available itkEMS [51] to segment rhesus brain. At the same time, some of the critical components utilized in these studies (e.g., registration tools) are not available, which does not permit the application of the developed pipelines to similar problems. It has also been shown that conventional tools originally developed for human MRI analysis can be applied to NHP imaging. McLaren et al. [15] used FSL tools for some of the processing steps to create the rhesus macaque atlas. Kochunov et al. used the FSL components to study fetal brain development in baboons [24,52]. Mietchen et al. [22] combined several tools including surface-based analysis suite by Van Essen et al. [53], FSL and custom filtering components for morphological brain analysis across multiple NHP species. These studies were made possible because of the public availability of the FSL tools. By developing and validating the presented processing pipeline, we show that 3D Slicer can also be used for the analysis of NHP MRIs. The added benefit of using 3D Slicer is the unrestricted ability to modify, customize and distribute its source code to meet any application-specific needs.

Our study has several limitations. Due to the specifics of the studied population, our atlas is representative only of the adolescent male vervet species. Our population was limited to 10 subjects, which may not be sufficient to represent anatomical variation. Our average intensity atlas is limited to T1 weighted (T1w) MRI. Although T2 weighted (T2w) sequences were included in the imaging protocol, we chose not to use these for atlas construction due to the significant discrepancy in slice thickness (0.5 mm slice thickness in T1w, 1 mm in T2w). The segmentation approach that we used can be further improved by incorporating shape constraints into the segmentation process, as proposed by Pohl et al. [48].

CONCLUSIONS

We presented an open source solution composed of the commonly available 3D Slicer tools to the automatic atlas-based segmentation of vervet brain MRIs. Our approach relies on computerized atlases of the average T1 signal intensity, and probabilistic atlases for the vervet brain MRI, which we developed as part of this work. The advantages of our pipeline over manual processing are reproducibility and reduced processing time. By carefully tuning the registration and segmentation steps, we produced high quality segmentations of challenging images, reducing the time needed for their processing, and improving the segmentation reproducibility. The flexibility of parameter tuning provided by the open source tools is critical for the analysis of non-conventional images. The animals used in our study derive from a colony of animals that are currently housed at Wake Forest University Primate Center. This colony is a NCRR resource that is part of a 9-generation, pedigreed and partially genotyped colony of vervet monkeys. As such, this animal resource will be exploited for multiple NHP studies and the presented methodology will be applicable to facilitate these studies in the near future.

We believe that the contributions of our work will assist in similar in vivo studies that rely on MRI for region-based morphological analysis of vervet brain. The digital vervet atlas we developed is now available as an open resource [32]. Our workflow is composed of free open source software tools available within 3D Slicer [27], which can be used without restrictions for either research or commercial purposes.

Table 1

Contrast to Noise Ratio (CNR) Between white Matter and Gray Matter in the T1w MRI Scans of the Studied Vervet Population. The CNR is Improved for Most of the Subjects in the Follow-Up Acquisition, and that there is Noticeable Variability in CNR Between the Individual Subjects

	GM/WM Contrast-to-Noise Ratio
Subject ID	1	2	3	4	5	6	7	8	9	10	mean (STD)
Baseline scan	9.81	13.42	16.89	9.52	5.24	13.81	9.63	7.23	10.94	8.66	10.51 (3.41)
Follow-up scan	13.73	15.06	14.61	14.87	13.6	11.72	13.5	9.08	14.07	17.18	13.74 (2.15)