A Magnetic Resonance Image Based Atlas of the Rabbit Brain for Automatic Parcellation

Rabbit brain has been used in several works for the analysis of neurodevelopment. However, there are not specific digital rabbit brain atlases that allow an automatic identification of brain regions, which is a crucial step for various neuroimage analyses, and, instead, manual delineation of areas of interest must be performed in order to evaluate a specific structure. For this reason, we propose an atlas of the rabbit brain based on magnetic resonance imaging, including both structural and diffusion weighted, that can be used for the automatic parcellation of the rabbit brain. Ten individual atlases, as well as an average template and probabilistic maps of the anatomical regions were built. In addition, an example of automatic segmentation based on this atlas is described.


Introduction
Animal models are essential for the understanding of brain and neurodevelopment. Several species have been used in neuroscience research, from primates to small animals such as rat, mouse and rabbit. The rabbit has been widely employed for modeling brain damage after perinatal injury in humans because it presents a human-like timing of perinatal brain white matter maturation [1]. Rabbit models of intrauterine inflammation [2], cerebral palsy [1,3] and intrauterine growth restriction [4] have been developed, demonstrating changes in neonatal neurobehavior and in brain structure [1][2][3][4][5][6].
Brain atlases have become an essential tool for the analysis of structural and functional differences in neuroimage, allowing volume and shape quantification of brain regions, for mapping functional activation and connectivity analysis. Over recent years, traditional 2D histological based atlases have been complemented by the generation of 3D digital atlases based on different image modalities, especially in magnetic resonance image (MRI). Although MRI-based atlases have less resolution than histological atlases, they present other advantages. Thus, 3D acquisition allows the volumetric reconstruction of brain regions, preserving the spatial relationship within the brain. Moreover, the digital format allows the application of image processing algorithms for quantification or automatic segmentation as well as comparisons between different subject acquisitions. Digital brain atlases have been developed for a number of species used in research, including mouse [7][8][9][10][11][12], rat [13,14], canary [15] or monkey [16][17][18][19].
However, to the best of our knowledge, there is no digital rabbit brain atlas available in the literature. MRI studies using the rabbit as a model have been based on manual delineation of the areas of interest. For this reason, we developed an MRI-based atlas for the New-Zealand rabbit brain, suitable for automatic segmentation. Delineation of regions was performed taking into account both T1-weighted and diffusion MRI, based on regions defined by histological atlases [20,21]. Nevertheless, some of the smaller regions described in these atlases, which cannot be identified radiologically, were not included in the template. The brain region delineation as well as the brain template and the probabilistic atlas is available on-line in www. medicinafetalbarcelona.org/rabbitbrainatlas, where the brain parcellation can be visualized and downloaded in order to be used for automatic segmentation.

Materials and Methods
In order to build the radiological rabbit brain atlas, T1 and diffusion MRI volumes of a set of 10 healthy adult rabbits were acquired and radiologically identifiable regions were manually delineated in these subjects. As a result, 10 individual brain atlases were obtained. Based on the 10 acquisitions, a brain template representing the average shape and intensity of T1-MRI brain volumes was built and a probabilistic atlas was developed, which defines at each point the probability of belonging to a specific region. Each of these steps are deeper described above.

Data
The atlas was constructed on a set of 10 healthy adult control New Zealand rabbits at 70 post-natal days (weight 2578+535 g, 40% male, 60% female). An additional healthy adult rabbit was used to test the performance of the region segmentation based on the atlas developed on the 10 experimental subjects. Animal experimentation of this study was approved by the Animal Experimental Ethics Committee of the University of Barcelona (permit number: 206/10-5440). Animal handling and all the procedures were performed following all applicable regulations and guidelines of the Animal Experimental Ethics Committee of the University of Barcelona. Included rabbits were obtained by Cesarean section at 30 days of gestation from New Zealand pregnant rabbits provided by a certified breeder. Rabbits were housed by a wet nurse rabbit until 30th postnatal day when they were weaned. Then, rabbits were housed in groups of three on a reversed 12/12 h light cycle with free access to water and standard chow. At 70th postnatal day, rabbits were anesthetized with ketamine 35 mg/kg and xylazine 5 mg/kg given intramuscularly and were sacrificed with an overdose of sodium pentobarbital (200 mg/kg) endovenous injection. Left and right common carotid arteries were cannulated and brains were perfused with phosphate-buffered saline (PBS) followed by 4% paraformaldehyde PBS. Finally, brains were dissected and fixed in 4% paraformaldehyde PBS at 4uC for 48 h.
The acquisition was performed on the excised and fixed brain using a 7 T animal MRI scanner (Bruker BioSpin MRI GMBH). High-resolution three-dimensional T1 weighted images were obtained by a Modified Driven Equilibrium Fourier Transform (MDEFT) 3D sequence with the following parameters: echo time (TE) = 3.5 ms, repetition time (TR) = 4000 ms, slice thickness = 0.7 mm with no interslice gap, 70 coronal slices and in-plane acquisition matrix of 188|188, resulting in a voxel dimension of 0:15|0:15|0:7 mm 3 .
For diffusion weighted images (DWI), Spin Echo DTI sequence was used to gain image quality, avoiding the artifacts associated to Echo Planar Imaging, but increasing the acquisition time [22]. Diffusion sensitizing gradients were applied along 126 directions with a b-value of 3000 s=mm 2 , and a reference (b~0) image was acquired. Other experimental parameters were: TE = 26 ms, TR = 250 ms, slice thickness = 0.7 mm with no interslice gap, 70 coronal slices and in-plane acquisition matrix of 40|40, with a voxel dimension of 0:7|0:7|0:7 mm 3 .

Image Processing
Previous to the manual delineation of the brain regions, image processing is required in order to take advantage of both T1 and diffusion MRI. The volumes acquired by both modalities were aligned, so T1 intensity and fiber orientation images can be jointly  visualized to perform the delineation. Since these modalities have different resolution, a multimodal registration algorithm was applied to align both images. Registration based on the optimization of mutual information [23] between T1 and the baseline volumes of the diffusion protocol was implemented. The affine transformation estimated by the registration algorithm was applied to the diffusion images, and afterwards the tensor image was estimated from the registered diffusion data set. The diffusion gradient direction is described with respect to the original image orientation. Consequently, changes in orientation due to the transformation applied to the diffusion images were also applied to the gradient direction [24].
In order to segment the brain from the background a mask was computed, by means of the Otsu threshold method [25]. Finally, the tensor at each voxel inside the mask was estimated using the least squares method described by [26].
Once the diffusion tensor image was computed, eigenanalysis was performed at each voxel. From eigenvalues, fractional anisotropy (FA) was computed and the first eigenvector was considered as the fiber direction [27]. Thus, the FA-color map, where color is related to fiber direction and intensity is weighted by FA was obtained.

Regions Definition
Taking as gold standard reference the histological rabbit brain atlas [20], manual delineation of brain regions was performed on T1-weighted images overlaid with FA-color maps. In addition, mouse and rat atlases [28][29][30][31][32] were used as second reference when structures were not described in rabbit atlas. Every brain structure was firstly delineated in the plane where was more clearly identifiable, and then corrected in the other two orthogonal planes. Although most of the regions were better identified in the coronal view, other planes were preferred for structures such as several cortical regions and the cerebellar vermis and hemispheres.
An example of the delineation of brain regions over representative slices of T1-weighted images is displayed in Figure 1. Furthermore, in the results section, the T1 intensity values and diffusion parameters characterizing each structure were compiled. 60 brain regions were defined, considering left and right structures separately when appropriate, which were classified into four groups: cortical regions, white matter (WM), deep gray matter (GM) and ''other regions'':  Note that region definition was based on radiological acquisitions, and therefore, finer regions requiring histological criteria to be identified are not included in the atlas. Without the aim of fully describe the delineated regions, below we include some guidelines taken into account to define the limits of certain structures, specially those structures that we have adapted from other species' brain atlases.   Regarding cortical regions, since no continuous cortical parcellation is available in the New Zealand's histological rabbit brain atlas [20], the delineation of most cortical areas, namely frontal, occipital, temporal, parietal, insular, piriform and entorhinal cortices, was performed based on mouse and/or rat atlases [28][29][30][31][32]. From anterior to posterior the cortical regions were labeled as follows: the medial portion of the cortex was defined as medial frontal cortex until the appearance of corpus callosum, after which was labeled as cingulate cortex. Frontal cortex region included the lateral parts of the cortex containing motor and sensory-motor areas [32]. The ventral part of the cortex was divided in olfactory, piriform and entorhinal cortices. Thus, following anterio-posterior direction, olfactory cortex was upper-limited by rhinal fissure. When rhinal fissure was not distinguishable, it became piriform cortex, which continued until the starting of the amygdala, where the beginning of entorhinal cortex was defined [28].
The delineation of WM regions was based on the work of Shek et al. [20]. Following anterior-posterior direction, we first found periventricular WM, which was considered as the WM surrounding lateral ventricles until the presence of the genus of the corpus callosum. Corona radiata, external and internal capsules were present also in the most anterior slices. When these structures, together with corpus callosum became not visible, subcortical WM is defined, until the end of WM bundles in the posterior part.
With regards to the GM regions, their delineation was based in the histological rabbit brain atlas [20] except for the amygdala, that was based in a rat atlas [28]. Namely, amygdala was identified  as the GM region surrounded by the posterior limit of the insular cortex and the anterior limit of entorhinal cortex. On the other hand, the anterior limit of the thalamic region coincides with the most anterior part of the fimbria of hippocampus. Thalamic region enclosed main thalamic and habenular nuclei, that would required histological analysis to be properly identified. The posterior limit of thalamic region was identified by the appearance of the superior colliculus. Hippocampus was easily identified as a multiple cortical layer structure in the coronal view and it included the hippocampal formation.
Finally, ''other regions'' category contained structures that did not fit in the previously define categories. This is the case of anatomical regions as cerebellar hemispheres, vermis, pons, medulla oblongata and septum and remainders of other brain regions as forebrain, basal forebrain, mesencephalon and diencephalon.

Delineation
The software used for delineation was ITK-SNAP [33]. It allows overlay of different images, with different transparency levels, and therefore delineation can be based on different image modalities. As aforementioned, both T1-weighted and diffusion magnetic resonance images were considered for a more accurate identification of the different structures composing the white and gray matter.
In order to simplify the delineation procedure, once the first image is delineated, its parcellation is propagated to the second subject by an elastic registration, so it can be taken as a starting point of the manual delineation of this volume as reported in [10]. This procedure is repeated iteratively to parcel the 10 subjects. At each step, all the previous delineations were considered, so a better starting point for the manual delineation is obtained.
Therefore, let be I 1 ,:::,I 10 the ten images to be parcelled and L 1 ,:::,L 10 the label maps corresponding to the parcellation of each of the subjects. Manual delineation of the the first volume I 1 resulted in a label map L 1 . Subsequently, every brain volume I n ,n~2,:::,10 was segmented based on the previous label maps L 1 ,:::,L n{1 , as follows: 1. The n{1 label maps previously obtained by manual delineation (L 1 n ,:::,L n{1 n ) were propagated to volume I n using an elastic registration algorithm. Thus, a set of n{1 label maps L L 1 n ,:::,L L n{1 n aligned to the volume I n were estimated. 2. A label map of subject n,L L n , was computed combining L L 1 n ,:::,L L n{1 n , assigning to each voxel x the most frequent label, that is:L L(x)~modefL L i n (x x)g, i~1,:::,n{1: 3.L L n is used as a starting point for the manual delineation of the brain regions of subject n, that results in L n .
This methodology resulted in a set of 10 individual atlases, that is, the region parcellation of the 10 brain volumes.

Average Template
A population template was built, describing the average shape and intensities of a normal healthy brain. The procedure followed to obtain this template was similar to the described in [34], first the average shape template is estimated iteratively, and afterwards the mean intensity model is computed: 1. Let be I n ,n~f1,:::,10g the ten volumes of healthy brains that were considered. 2. The most normal volume, I min disp in the data set was chosen to initialize the iterative algorithm. It is defined as the volume requiring the minimum transformation to match all the other volumes in the dataset. The elastic transformation matching every pair of volumes was estimated by means of a block matching registration algorithm [35], resulting in a displacement vector field for each pair of images. For each of these transformations the mean displacement was computed. Finally, the image minimizing the mean displacement was used to initialize the iterative algorithm followed to determine the average shape template. 3. Once the minimum displacement image was identified, the mean shape, T shape was computed by means of an iterative procedure. Let be T 0~Imin disp the initial estimation. It was registered against all the images I n in the dataset, and the average transformation T 0 was computed. This transformation was applied to the current template to obtain the template for the next iteration T iz1~T i (T i ). This procedure was repeated until the average transformation was smaller than a given threshold. In practice, convergence was achieved in few iterations. The mean shape template was obtained as T shape~T i Ã (T i Ã ), where i Ã is the iteration in which convergence is reached.
4. Taking into account the mean shape, the average intensity volume was computed. That is, all the volumes were registered to the mean shape and the average intensity value at each voxel was computed. Voxels whose intensity was above two standard deviation of the mean value were excluded to avoid the effect of noise or misregistration in the template.

Probabilistic Atlas
A probabilistic atlas was built over the template based on the 10 individual atlases, describing at any location the probability to belong to any of the regions. Dice coefficient between the manually delineated brain regions and the brain regions identified by the automatic atlas-based segmentation, and between the manual delineations performed by two different observers. doi:10.1371/journal.pone.0067418.t002 The label maps L n , n~1,::,10 that had been delineated over each volume I n ,n~1,:::,10 were propagated to the average template, resulting in a set of ten label maps L avg n . At each voxel x x of the average template, the probability of belonging to a given region k was estimated as where D : D denotes the cardinal of the set, and N is the number of volumes considered to build the atlas, that is, N~10. Thus, we obtained a set of probabilistic maps, one for each anatomical region delineated in the atlas. The use of this probabilistic approach is more robust against volume partial effect, since voxels in the edge between two regions (let be R i and R j ) will have a certain probability P i to belong to R i and a probability P j to belong to R j , which is especially useful for the automatic parcellation. Also a label map can be estimated on the template assigning to each voxel the label of the most probable region.

Automatic Parcellation
The atlas can be used for automatic brain parcellation based on registration. Let be I a new brain volume, segmentation is obtained by registering the template T against it, assessing in that way the elastic transformation T : T?I. This transformation can be estimated by any of the software available for image registration. Applying this transformation to the region probability maps, the probability of a voxel in the image I to belong to each of the regions is computed. Finally, each voxel is assigned to the region of maximum probability. It is also feasible to apply the transformation to the label map defined over the average template, obtaining in such way the label map in the new brain, although it could be less accurate than the probabilistic approach.
On the other hand, the accuracy of the segmentation relies on the performance of the registration algorithm. To test the approach, a multiresolution block-matching algorithm was implemented to perform registration, based on the correlation coefficient between T1-images [35]. The performance of segmentation is tested in the additional subject that was not included in the atlas building, and evaluated both qualitatively, by visual inspection, and quantitatively, by comparing with manual delineations performed by two different observers. Namely, Dice coefficients and confusion matrix [36] were used to measure the similarity between manual and automatic segmentation. The overlapping between two different parcellations was estimated by the Dice coefficient for each region i: were D : D is the cardinal of the set, X i are the points that were labeled as belonging to region i by the first parcellation being compared and Y i the points assigned to region i by the second parcellation. This index is computed to measure the similarity between automatic and manual delineation, as well as between both manual delineations. High and similar values of Dice coefficient in both cases will show the reliability of the automatic segmentation, meaning that differences are comparable to the interobserver variability.
In addition, a measure of the global matching was estimated as: Besides, confusion matrix, measuring the percentage of voxels belonging to region i that have been labeled as region j, was built.
Obviously, automatic segmentation can be performed with other registration algorithms [37], and the accuracy of the result will rely on the performance of the registration algorithm.

Anatomical Regions
As previously described, a set of 60 brain regions was defined for each volume. Each region was assigned to one of the following areas: cortical, white matter, deep gray matter, and ''other regions''. Illustrative views of the 3D reconstruction of the regions included in the four main areas are depicted in Figure 2. Figure 3 displays representative slices of the T1-images, where the corresponding regions are overlapped for each major area. The properties of the structures, such as T1 intensity and diffusion parameters, are compiled in Table 1, as described below.

Individual Atlases
Ten individual atlases were developed as described in the Methods section. In Figure 4, different slices of the individual atlases are shown. Each row corresponds to a different subject and each column represents equivalent slices for each of the subjects, containing similar structures.
Different regions are characterized in Table 1, where mean values of regional volume (corrected for total brain volume) and relative regional T1 intensity (normalized by average T1 intensity in the whole brain) of brain regions are displayed, together with the mean values of regional fractional anisotropy (FA) and regional mean diffusivity (MD) normalized by average MD value in the whole brain. Thus, relative values will be higher than 1 if they are higher than the mean value in the whole brain, and lower than 1 if they are lower than the mean value.
Note that separate left and right sides of most bilateral structures had been taken as different regions in the atlas, obtaining 60 regions. However, for the sake of simplicity, in Table 1, left and right sides were considered altogether, resulting in 35 different regions.

Template and Probabilistic Maps
Some slices of the template volume are shown in Figure 5. The probability maps of some of the regions are shown in Figure 6, where it can be noted that the contours of regions are fuzzy, since voxels in these areas may belong to neighbor regions.

Automatic Segmentation
In order to test the accuracy of segmentation based on the atlas, the brain volume not included in the atlas building was automatically segmented, and compared with the manual delineation of this volume. Segmentation performance was both qualitative and quantitatively evaluated. First, visual inspection of the resulting segmentation confirms appropriate segmentation, as can be viewed in Figure 7, where some slices of the T1-MRI and the overlapped contours of the automatically segmented regions are displayed. It can be observed that the different structures were correctly identified, even in areas where the tissue is broken.
Secondly, objective measures also confirm the similarity between the manual and automatic segmentation: the index for global matching between the automatic segmentation and each of the two manual delineations were 0.9187 and 0.8690; and the Dice coefficient computed between both manual delineation was 0.8779. That is, globally, the accuracy of the automatic parcellation is similar to the accuracy of the manual delineations.
The accuracy of segmentation for each individual region is compiled in Table 2. Note that right and left areas of the same structure are considered as an only region. The three columns in table correspond to: similarity between automatically identified regions and the first manual delineation; similarity between automatically identified regions and the second manual delineation; and similarity between both manual delineations.
Finally, the confusion matrix is shown in Figure 8. The value at each point (i,j) in the matrix is the percentage of voxels belonging to region i in the manual delineation that have been labeled as region j by the automatic segmentation. That is, brighter points corresponds to higher number of points belonging to region i labeled as region j. In case of perfect matching, diagonal values would be one (white) and the others point would be zeros. It can be viewed that the resulting confusion matrix for automatic segmentation is close to be diagonal.

Anatomical Areas Definition
Having an MRI-based rabbit brain atlas to allow automatic segmentation is of great interest since it opens a wide window for neuroimage based analysis as, for instance, connectivity studies. In this regard, manual delineation was performed using both T1-MRI and diffusion MRI data. This multimodal approach allows a more accurate identification of specific structures such as WM tracts. However, the spatial resolution of both types of images limits the delineation of different anatomical structures. For this reason, all structures that could not be delineated were distributed into major divisions of the central nervous system, such as those described in ''other regions'' category.
Note that delineation was performed over images of postmortem fixed and excised brains. It must be taken into account that there are morphometric differences between in vivo and in vitro brains [11]. For this reason, special care must be taken if the atlas is applied to segment images of in vivo brains, being necessary appropriate registration algorithms to remove the postfixation distortion.

Individual Atlases
Ten individual brain atlases were built in order to avoid the bias due to the choice of a single subject. The low variability among size and intensity values in the ten subjects supports that parcellation of brain regions was highly reproducible.
The delineation scheme here used has been already reported in [10], and simplifies the tedious task of manual delineation.
Although it could be argued that a previous automated delineation step can bias the observer, we countered this potential drawback by manual correction verifying the compliance of the delineations with rigorous criteria as assessed by an expert in neuroanatomy.

Template and Probabilistic Maps
The use of an average model of shape and intensity allows to have a reference template which represents the normal shape and intensity distribution of the rabbit brain, avoiding the inter-subject variability. Probability maps are used in order to deal with the partial volume effect. While individual atlases assign single values to each voxel to identify the region to which the voxel belongs, the probabilistic maps give a probability. This approach may allow a higher accuracy in the definition of region when a new sample image is registered to the template. The number of subjects required to obtain a probabilistic map is not clearly defined. Previous studies suggest that the use of 10 subjects as performed in this study allows to build representative probabilistic atlas [10].

Automatic Parcellation
In this paper, we have proposed an automatic segmentation method based on the maximization of the region probability at each voxel. To match the template to the data image, a multiresolution block-matching algorithm based on the correlation coefficient between the intensity levels was used. The use of this algorithm allows robust global matching avoiding local minima. However, the atlas here reported could potentially be used with other registration algorithms, such as the implemented in available image processing software.
Quantitative evaluation showed that the differences between the regions automatically and manually identified were comparable to the differences due to the interobserver variability (Table 2), which supports that the atlas can be used for automatic brain parcellation in studies using the rabbit brain. All the regions could be automatically identified by means of registration against the proposed atlas with accuracy values similar to the interobserver differences. It can be noticed that similarity values were higher in bigger regions that in smaller nucleus. In these smaller areas, subtle differences in the contours of the regions have more influence in the final measure of the Dice coefficient, since they represent a higher percentage of all the voxels belonging to the region. For this reason, lower similarity values in these regions were present in the comparison between automatic and manual regions as well as in the comparison between manual delineations. The only region where there was a significant improvement when delineation was performed manually was the olfactory bulb. This fact could be related to the high variability of this structure in our data-set, due to the brain extraction and fixation procedure.

Conclusions
Atlases have become fundamental in neuroimage, since they are required to identify brain structures in a coherent and objective way in different subjects. Moreover, the use of digital atlases allows automatic segmentation of such structures, avoiding the necessity for manual delineation to perform regional analyses. In this paper, we contribute to solve the lack of digital atlases of the rabbit brain by developing an MRI-based atlas of the New-Zealand rabbit available on line. First, a set of anatomical regions that constitute the rabbit brain have been defined based on the literature. These regions have been identified in a set of ten individuals, showing the reproducibility of the anatomical parcellation in different subjects.
One of the main applications of the anatomical atlas here described is to be used for automatic segmentation. An average template and a probabilistic atlas have been developed from the individual atlases in order to provide a subject-independent reference of brain parcellation and a model of normality for the brain. Moreover, the template and the probabilistic atlases are useful for the development of automatic segmentation algorithms. The ability of the atlas to be used for automatic segmentation has been tested, and the quantitative comparison with manual delineation has shown that similar results are obtained.
Therefore, the atlas here presented will be a useful tool for studies using the rabbit as a model of brain disease.