Myelo- and cytoarchitectonic microstructural and functional human cortical atlases reconstructed in common MRI space

The parcellation of the brain's cortical surface into anatomically and/or functionally distinct areas is a topic of ongoing investigation and interest. We provide digital versions of six classical human brain atlases in common MRI space. The cortical atlases represent a range of modalities, including cyto- and myeloarchitecture (Campbell, Smith, Brodmann and Von Economo), myelogenesis (Flechsig), and mappings of symptomatic information in relation to the spatial location of brain lesions (Kleist). Digital reconstructions of these important cortical atlases widen the range of modalities for which cortex-wide imaging atlases are currently available and offer the opportunity to compare and combine microstructural and lesion-based functional atlases with in-vivo imaging-based atlases.


Introduction
Parcellation of the cortex into distinct and discrete subareas is a topic of ongoing interest ( Eickhoff et al., 2018 ). Pioneering early-day neuroanatomists such as Campbell, Smith, Brodmann, Von Economo, Flechsig and Kleist defined regional boundaries based on heterogeneity of cortical microstructure ( Brodmann, 1909 ;Campbell, 1905 ;Flechsig, 1920 ;Smith, 1907 ;von Economo and Koskinas, 1925 ) and functional segregation ( Kleist, 1934 ). These cortical mappings are based on extensive studies of post-mortem brain tissue using different staining methods, revealing local variation in cyto-and myeloarchitecture. Detailed documentation of the spatial location of cortical lesions and resulting clinical symptoms formed the basis of early functional mappings of the human cerebral cortex ( Kleist, 1934 ). Since their introduction over a century ago, region descriptions from such neuroanatomical brain atlases are widely used in neuroscientific research ( Galaburda et al., 1978 ;Shafiei et al., 2020 ).
We provide digital reconstructions of the Campbell ( Campbell, 1905 ), Smith ( Smith, 1907 ), Brodmann ( Brodmann, 1909 ), Flechsig ( Flechsig, 1920 ), Von Economo -cortical type ( von Economo and Koskinas, 1925 ) and Kleist ( Kleist, 1934 ) Table 1 Demographic overview of classic brain mapping studies. Demographic overview of the six digitized atlases, describing the year of publication, documented number of areas in the original and digital atlases, modalities upon which the atlases were based, documented number of hemispheres on which the atlas was based, and the age of the examined subjects, A ( Brodmann, 1909 ), B ( Flechsig, 1920

Atlas descriptions
This section provides an overview of the original atlases and their digitized counterparts, describing per atlas the modality examined and the number of included regions, together with relevant technical details of the methods through which the original atlas was constructed. A demographic overview of the number of areas, the modalities used for defining regional boundaries and the number of hemispheres used is given in Table 1 . Tables S1-S5 describe for each atlas the label names of the regions as included in the digitized atlases and their corresponding area nomenclature from the original atlas description. All brain atlases describe an average hemisphere, with no further detail on left-right asymmetry, following the documentation of the original atlases.

Campbell 1905
Campbell studied the histological localization of cerebral function based on cyto-and myeloarchitecture and constructed a parcellation of the cortex using regional variation of the size, laminar arrangement and number of fibers and neurons as criteria for regional boundaries. The 1905 version of the Campbell atlas describes 17 distinct cortical areas, with notes on possible functions based on experimental, developmental, clinical and comparative information ( Campbell, 1905 ). Regions were named according to their function and/or location such as visual-sensory area and temporal area ( Campbell, 1905 ).

Smith 1907
The 1907 Smith atlas describes 43 distinct cortical areas ( Smith, 1907 ). Regional boundaries were defined based on local variation in cyto-and myeloarchitecture and cortical thickness. Main sulci and gyri were described, ignoring minor gyrifications more variable across the population ( Smith, 1907 ). Brain regions were named after the main sulci and gyri within which a region was located, including a suffix when further subdivision within a sulcus/gyrus was observed (such as area superior frontal A).

Brodmann 1909
Brodmann's histological studies resulted in a widely used cortical parcellation scheme of 47 cortical areas ( Brodmann, 1909 ). Areal boundaries were defined based on regional variation in cell size and composition, cortical thickness and composition of cortical layering. Each brain region was labeled with a unique number, with numbering starting at the central sulcus. We merged from these 47 regions the following sets of small and/or oblong regions to ensure cross-population delineation consistency of the digitized atlas: Brodmann area BA1 and BA3, two narrow and oblong areas located on the posterior wall of the central gyrus (merged into BA1_3); Brodmann areas 26, 29 and 30, three small and narrow areas in the posterior cingulate (merged into BA26_29_30); Areas 41, 42 and 52 (Heschl's gyri) located in the bank of the dorsal surface of the temporal lobe (merged into BA41_42_52); Area 35 and 36 representing the entorhinal and perirhinal cortex (merged into BA35_36). Area 48 (subdivision of the hippocampus, defined as a subcortical area) was excluded from the digital cortical atlas. Please note that as Brodmann reported areas 12,14,15,49,50 and 51 to be present in other species, but not found to have an equivalent in the human brain, these regions were not included in the Brodmann human cortical atlas ( Brodmann, 1909 ).

Flechsig 1920
Flechsig studied myelogenesis in developing brains from fetuses to adults ( Flechsig, 1920 ). Based on the temporal profile of myelination, Flechsig parcellated the cortex into 54 different areas, including nine subdivisions of cortical regions indicated with a letter as suffix (e.g. area 2, area 2b). A numbering system was used to specify each region, with lower numbers representing cortical areas being myelinated early in development and higher numbers describing areas myelinated at a later time point. For digital reconstruction a few modifications were made relative to the original atlas: Area 1 was described to cover the lower part of the central sulcus ( Flechsig, 1920 ), but we were unable to find further detail on the exact boundaries of this area in the original documentation. We therefore took the bottom of the central sulcus as the boundary between areas 3 and 4 (post-and pre-central regions, respectively). Further subdivisions of area 3 (3b), 4 (4b and 4c), 5 (5a and 5 ′ ), 14 (14b) and 18 (18a) ( Flechsig, 1920 ) were not accounted for in terms of separate labels in the digital atlas as we could not find details on the exact spatial location of these boundaries. These areas were thus merged into their larger main area.

Von Economo 1925 -cortical type
The Von Economo and Koskinas atlas describes five broad fundamental structural types and one intermediate, or transitional, type of cortex based on histological observations of laminar organizational patterns ( von Economo and Koskinas, 1925 ); their cytoarchitectural atlas has been digitized in a previous study . Cortical areas were subdivided into five cortical classes (cortical types) ( von Economo and Koskinas, 1925 ): Cortical type 1 (agranular cortex) was observed to lack granule cells and to thus have no recognizable layer II and IV. Cortical type 2 (frontal), 3 (parietal) and 4 (polar) were classified as homotypic cortex and named after their most prominent spatial location. Homotypic cortex was characterized by a conserved six-layered organization ( von Economo and Koskinas, 1925 ). Cortical type 5 (granulous cortex) was described to contain large numbers of granule cells, even in layers III, V and/or VI. Transitional cortical type 1/2 was described as areas with similar cortical thickness and cell morphology to type 1, but to contain granular layers ( von Economo and Koskinas, 1925 ). We segmented these cortical types, with non-contiguous surfaces of each cortical type implemented as separate labels.
2.1.6. Kleist 1934 The Kleist atlas describes anatomical and behavioral data of patients suffering from traumatic brain injuries and nontraumatic focal brain lesions ( Kleist, 1934 ). Through the loss-of-function resulting from these lesions, Kleist deduced localized cerebral functions ( Kleist, 1934 ;Neumärker and Bartsch, 2003 ). Kleist defined the exact spatial location of the brain lesions post-mortem and linked them to the previously documented loss-of-function of the patient. For his atlas, Kleist used Brodmann's definition of areas as a template to superimpose his functional mapping on ( Kleist, 1934 ). Similarly, we used the digital Brodmann atlas ( Section 2.1.3 ) as a basis for the digital Kleist atlas. Kleist used Brodmann areas as anatomical landmarks and included further functional subdivisions within certain areas where possible (for example, Brodmann area 6, 17, 18, 19, 22, 39, 40 and 44). Brodmann areas 28 and 34 were assigned a common function and were thus merged in the digital atlas. Cingulate regions (Brodmann area 23,24,26,29,30,31 and 32) were also allocated with a common function and were combined into a single label in the digital atlas. Reported specification of functional localization of face, trunk and extremities within area 4 and 6 (post-and precentral) as reported by Kleist was not further included as subareas in the digital atlas, as we were not able to find detailed documentation on the exact boundaries to delineate them separately onto the digital template brain. In addition, we could not find specific documentation on the functional allocation of Brodmann areas 5,13,16,25,27 and 38. These regions were therefore labeled with 'no assigned function' in the digital atlas. This resulted in a total of 43 brain regions (with six additional labels representing the six Brodmann areas without functional allocation) included in the digital atlas. A list of the assigned functions per region as defined by Kleist (translated from German to English) is provided in Table S6 ( Kleist, 1934 ).

Atlas digitization
This following section provides an overview of the process of manual segmentation used to digitally reconstruct the original atlases ( Fig. 1 ), including 4 main steps:

Step 1. Standardized space
The Colin27 template brain of the McConnell Brain Imaging Centre (BIC) of the Montreal Neurological Institute of McGill University ( Aubert-Broche et al., 2006 ;Holmes et al., 1998 ) was used as a segmentation template. The Colin27 average brain was built from 27 T1-weighted scans of the same subject, with the purpose of facilitating detailed segmentation of the cortex ( Evans et al., 2012 ;Holmes et al., 1998 ). The Colin27 template was processed using the standard FreeSurfer pipeline ( http://surfer.nmr.mgh.harvard.edu/ ; version 5.3.0), resulting in a high-resolution surface reconstruction of both hemispheres ( Fischl, 2012 ).

Step 2. Segmentation of brain regions
Described regions were segmented on the cortical surface of the left and right hemisphere of the Colin27 brain for each atlas, based on the original illustrations and descriptions of each region by the original authors. The digital atlases each describe an average hemisphere, with no further detail on left-right asymmetry, following the documentation of the original atlases. To enable a 2D-to-3D translation of the schematized cortical representation of the original atlas, the major gyri and sulci were matched to those of Colin27, forming the starting point of the manual digitization process. Then, each individual region was positioned relative to these major landmarks and manually segmented using the tksurfer label segmentation and editing tool included in the FreeSurfer package Fischl et al., 1999 ). Each region was delineated by plotting individual landmark points onto the cortical surface using tksurfer 's drawing tools. In those locations where a region boundary was situated on the lateral surface of the cortex, the . pial surface view was used, while in those instances where the region boundary continued into the depth of a sulcus, the view was switched to the . white or . inflated surface (overlaid with curvature data showing sulcal depths and gyral peaks for orientation purposes) ( Fig. 1 ). In cases where region boundaries did not coincide with major landmarks, the boundary was plotted to provide a gyral subdivision proportional to that shown in the original atlas illustration. Following this process, for each region, landmark points were plotted along the full region boundary and connected. Finally, the closed boundary was filled with the paint label tool to create a color coded surface for export as a . label file.

Step 3. Construction of annotation files
Labels were visually inspected and corrected for gaps and/or overlap between labels across all regions. The position of the manually segmented labels was compared to the original descriptions and illustrations of the respective atlas and adjustments were made where necessary. Label segmentations were performed by RP in close collaboration with LHS to ensure good segmentation quality. A .ctab file containing a list of RGB color codes for all individual labels was created. Next, in each hemisphere, mris_label2annot combined the individual .label files together with the .ctab color codes, creating a .annot annotation file listing the vertex label assignments for all regions and their respective label color.

Step 4. Atlas formation
The annotation file was used as input for the mris_ca_train function to generate an atlas ( .gcs ) file per hemisphere, generated by combining the spherical surface registration from Colin27 to the FreeSurfer fsaverage subject (to which all individuals processed with the standard FreeSurfer pipeline are registered) with the manually segmented labels included in the .annot file. Atlas label consistency was visually evaluated on datasets from the Q3 release ( n = 500, T1-weighted, 0.7 mm isotropic resolution) of the Human Connectome Project (HCP) ( Glasser et al., 2013 ;Van Essen et al., 2012 ), using the mris_ca_label function to create individual annotation files for each of the atlases. A random subset of 20 subjects was examined in more detail to report on the general applicability of the atlases to external MRI data and to make small corrections to the final versions of the atlases (using this information again in step 3 and 4 of the digitization process). Example code for application of the surface-based atlases to datasets of other FreeSurfer-processed subjects is provided in the Supplemental Methods.

Transformation to fsaverage GIFTI format
To facilitate incorporation in surface-based processing steps outside of the FreeSurfer environment, for each of the six atlases, the two hemisphere-specific surface atlas ( .gcs ) files were applied to FreeSurfer's fsaverage subject, generating two single-hemisphere annotation ( .annot ) files per atlas. Next, each of these .annot files were converted into GIFTI format label files using the mris_convert function in FreeSurfer, resulting in a set of twelve ' l/rh.fsaverage.atlas.label.gii' files.

Transformation to volumetric coordinate system
The six digitized surface-based atlases were transformed to volumebased atlases for the Colin27 template and the MNI ICBM 2009a Nonlinear Symmetric (hereafter referred to as ICBM152) stereotaxic template for incorporation in voxel-wise analyses ( Fonov et al., 2009 ). Similar to the analysis of the Colin27 template, the ICBM152 template was processed using the standard FreeSurfer processing pipeline ( http://surfer.nmr.mgh.harvard.edu/ ; version 5.3.0). Annotation files ( .annot ) were generated for Colin27 and ICBM152 using the surfacebased atlas files ( .gcs ) for both the left and right hemisphere. Resulting annotation files were then used as input for the mri_aparc2aseg function (mapping the automatic surface parcellation to the aseg automatic segmentation volume) to convert the surface-based atlases to volume-based atlases. The voxel-label assignment was visually inspected to ensure consistent label placement across atlas modalities. Fig. 3 shows the Colin27 volume-based version of the digitized Campbell atlas (see Supplemental In this overview of the atlas digitization process, the left hemisphere and the Campbell atlas is taken as an example ( Campbell, 1905 ). A) The left panel represents the original documentation of the classic atlases, including an image of the brain map. The right panel shows the lateral and medial view of the Colin27 reference brain, used as a segmentation template. B) The tksurfer tool, included in the FreeSurfer software, was used for manual segmentation of the brain regions on the Colin27 surface (right). The first region to be delineated is highlighted in red in the original drawings on the left. To delineate a brain region, individual landmark points were plotted on the cortical surface (white dots) after which these landmark points were connected (red line). The cortical area was filled and saved as a label ( .label ) file. C) The same procedure iterates for the next region (left panel, in yellow). All brain regions were segmented this way. Across all regions, labels were then merged into an annotation ( .annot ) file. D) . During the 3D delineation process, the surface visualization was alternated between the 'normal' pial surface (top), the inflated surface (middle), and the white matter surface (bottom). The annotation file was translated to an atlas ( .gcs ) file, forming a cortex-wide brain atlas. For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Figures S10 and S11 for the other atlases). The volume-based atlases are provided as nii.gz files together with an atlas-specific lookup table.

Atlas consistency
Atlas consistency was evaluated using three analyses. We first compared the spatial consistency of region labels across the set of six atlases. Second, we compared the region labels with the commonly used surface-based ex vivo Brodmann areas included in FreeSurfer ( Fischl et al., 2007 ) and the benchmark probabilistic volume-based region labels of the Julich-Brain cytoarchitectural atlas ( Amunts et al., 2020 ). Third, we overlaid each atlas with NeuroSynth functional annotations ( Yarkoni et al., 2011 ) and examined the atlas region -functional term overlap.

Spatial consistency across atlases
Spatial correspondence of region labels and boundaries across the set of digitized surface atlases was evaluated by computing atlas overlap on three levels: 1) global correspondence across all six atlases; 2) vertex-wise spatial consistency of region assignment; 3) vertex-wise consistency of region boundary placement.
Atlas overlap was computed between all atlas pairs, with the overall consistency of atlas A with atlas B computed as the shared surface area of each region across the atlases. For a region i in atlas A (region A i ), the consistency with region j in atlas B (region B j ) was the proportion of shared surface area: For each region in atlas A, the most consistent region in atlas B was selected, and cross-atlas consistency was computed as the average region consistency across all regions of atlas A , resulting in a global consistency score per atlas pair ( Fig. 4 A, Supplemental Table S7, Supplemental Figure S8). This asymmetric consistency metric is a variant of the intersection over union metric and has the benefit that it is also suitable for comparing atlases in which cortical regions differ in size between atlases. For example, this metric is in line with the intuition that in the event that atlas A is a subdivision of atlas B (e.g. the case for Kleist and Brodmann), the region consistency of A to B will be close to 1 and B to A will be much smaller than 1.
To determine where across the cortex the six atlases showed the greatest consistency, vertex-wise spatial consistency of region assignment across all atlases was computed. For each atlas, for each region, the region consistency with the other five atlases was computed (as above) and averaged,yielding six atlas consistency maps, depicting the extent to which vertex-region assignment of a given atlas ( A ) was consistent with each of the other five atlases ( B-F ). Next, for each vertex in the surface file, region assignment consistency was averaged across atlases, resulting in a hemisphere-specific map illustrating areas of high vs low consistency across atlases ( Fig. 4 B; Supplemental Figure S6).
Vertex-wise consistency of boundary placement across atlases was computed by counting the number of times a vertex was involved in a regional boundary as indicated by having a neighbor vertex with a different label. Boundaries involving regions labeled as unknown (segmenting out non-cortical structures on the medial wall of the surface reconstruction) were excluded from the boundary consistency measure. This resulted in a map visually highlighting where on the cortical surface region boundaries tend to coincide across all six atlases ( Fig. 4 C, Supplemental Figure S7) as well as between all atlas pairs (Supplemental Figure S9).

Spatial overlap with probabilistic cytoarchitectonic atlas labels
We further examined the spatial similarity in label placement between the presented atlases and cytoarchitectonic region labels. We computed the vertex-wise overlap of eleven surface-based probabilistic ex vivo Brodmann areas ( Fischl et al., 2007 ) with matching regions of the digitized Brodmann surface atlas (see Supplemental Table S8). After application of the ex vivo labels to Colin27 (using FreeSurfer's recon-all -s subjid -ba-label s), for each of the ex vivo labels, the list of vertices belonging to the ex vivo label and its matched digitized label was extracted, and the proportion of shared vertices was computed, quantifying the extent to which both labels were co-localized on the cortical surface reconstruction.
We further compared the probabilistic region labels of the Julich-Brain cytoarchitectonic atlas ( Amunts et al., 2020 ) with the Colin27 MNI volume-based versions of the digitized atlases. Colin27 hemispherespecific probabilistic Julich-Brain region labels (110 per hemisphere) were downloaded from the Julich-Brain atlas viewer (at https://julichbrain.fz-juelich.de/ ), specifying for each region the voxel-wise probability that it is located at a given position on the cortex. The level of voxelwise region overlap with all Julich-Brain labels (probability > 0.1) was computed for each of the digitized atlases, calculating the ratio of voxel overlap per region (i.e., the number of overlapping voxels divided by the total number of voxels within a region), together with the mean Julich-Brain probability score of those overlapping voxels. The most overlapping Julich-Brain regions (with a minimal overlap of 10%) were ranked, providing an overview of the level of spatial similarity between the digitized atlas regions and the probabilistic Julich-Brain atlas.

Spatial overlap with NeuroSynth functional annotations
We quantified the spatial overlap of each atlas' respective regions with functional annotations (i.e. meta-analysis images of manuscript text mined terms and study ROIs) to map the cortical involvement in a large variety of cognitive tasks and brain states as documented in the comprehensive NeuroSynth database ( http://www.neurosynth.org ) ( Yarkoni et al., 2011 ). All brain-related functional annotations ( N = 942, representing voxel-wise z-scores quantifying the extent to which a voxel is associated with a corresponding functional term; non-significant voxels were set to zero, FDR corrected, q < 0.05) were downloaded from the NeuroSynth online database and registered to the ICBM152 template using FLIRT ( Jenkinson et al., 2002 ;Jenkinson and Smith, 2001 ). For each cortical region of the six digitized atlases, for each NeuroSynth term, the ratio of overlapping voxels (i.e., the number of overlapping voxels divided by the total number of voxels within a region) was computed, reflecting the extent to which that region in the volume-based ICBM152 atlas was involved in the NeuroSynth term. Given that the NeuroSynth maps are bihemispheric, contralateral regions of the digitized atlases were combined into a bihemispheric region of interest for analysis. For each atlas region, the top 10 most matching terms were selected, together describing a putative functional profile of that region (Supplemental Table S10).
We additionally investigated the spatial overlap of the region functions described by the Kleist (1934 ) lesion-based functional brain map and NeuroSynth terms. For each region, the conceptually best-matching NeuroSynth term was selected, after which from the voxels belonging to that region, the subset with a significant functional activation related to the NeuroSynth term was extracted. Finally, the mean z -score across the extracted voxels was calculated to represent the extent to which there was a significant spatial correspondence between the region label and its matched term (Supplemental Table S11). Fig. 2 shows illustrations of the original atlases side by side with the corresponding digital atlases superimposed on the Colin27 brain template. The digital reconstruction of the Campbell atlas ( Fig. 2 A) describes 17 cortical areas, with each color representing a different area. For the digital reconstructions of the atlases of Smith ( Fig. 2 B; 43 regions), Brodmann ( Fig. 2 C; 39 regions) and Flechsig ( Fig. 2 D; 46 regions) respectively, areas are color-coded by lobe, with different shades within a lobe. The six cortical types included in the Von Economo cortical type atlas ( Fig. 2 E) are segmented in 15 label files (Table S5), with labels of cortical type 1 through 5 presented in red, yellow, green, light blue and dark blue respectively, intermediate type 1-2 is depicted in orange. For analyses, anatomically separate labels within each cortical type may be merged. The color coding by lobe of the digital Kleist atlas ( Fig. 2 F; 49 labels) is similar to Brodmann's digital atlas, with cortical areas that had no allocated function in shades of gray. Table S6 provides an English translation of Kleist's functional allocations. Tables S1-S5 describe a list of all regions included in the other cortical atlases, Figures S1-5 show surface views of the atlases applied to five randomly selected HCP subjects.

Download packages
Packages for all six atlases are included as supplemental material. All atlases, as well as any future versions are also available for download at www.connectomelab.nl/ (downloads). Each atlas package includes: (1) a README file, describing the version, date, agreement of use and Creative Common License. The README file includes a table describing the regions in the atlas and the corresponding digital labels. For the Kleist atlas, regional functional allocation is included in original German and translated English. (2) A folder containing the individual manually created label files. (3) Surface-based atlas files ( lh.atlas.gcs and rh.atlas.gcs ) for the left and right hemisphere. (4) Corresponding color table files, listing individual labels with respective RGB color coding ( lh.colortable.txt and rh.colortable.txt ). (5) GIFTI label files in fsaverage surface space. (6) Volume-based atlas files for the Colin27 brain ( ATLAS.nii.gz ) and the ICBM152 template ( ATLAS_ICBM152.nii.gz ) with corresponding (7) lookup table ( ATLASColorLUT.txt ) containing ROI numbering and color coding.

Spatial consistency across atlases
Average global consistency between atlas pairs was 0.65 (std 0.12), with the lowest consistency observed at 0.49 (Von Economo cortical type -Kleist) and the highest 0.97 -as expected-observed between the Brodmann ( Brodmann, 1909 ) and Kleist atlas ( Kleist, 1934 ), with the  Kleist and their digital counterparts. The 3D digital atlases are shown on the cortical surface reconstruction of the Colin27 template brain on which they were delineated. ( Brodmann, 1909 ;Campbell, 1905 ;Flechsig, 1920 ;Kleist, 1934 ;Smith, 1907 ;von Economo and Koskinas, 1925 ).  ( Campbell, 1905 ) ( Campbell.nii.gz in the atlas download package), as visualized using Freeview ( Fischl, 2012 ). A shows the location of the slices visualized in B (sagittal), C (coronal) and D (transverse). The individual regions are labeled using the same color table as in the surface atlas (Supplemental Figures S10 and S11 show all volume-based atlases on Colin27 and ICBM152 respectively). majority of Kleist regions being subdivisions of the Brodmann atlas ( Fig. 4 A, Supplemental Table S7).
Quantifying the degree of vertex-wise consistency in region assignment across atlases showed that the areas with the highest consistency are located in primary visual and motor cortex, and to a somewhat lesser extent in the medial temporal and orbitofrontal cortex, while areas with lower consistency are located in superior frontal, superior parietal and cingulate cortex ( Fig. 4 B, right hemisphere and inflated view in Supplemental Figure S6).
Vertex-wise consistency of region boundary placement shows a similar pattern, with major boundaries delineating distinct cortical types (e.g. between V1 and V2, primary motor and sensory cortex; or delineating cingulate or insular cortex) highly consistent across atlases, while boundaries in e.g. prefrontal and superior parietal cortex are much more variable ( Fig. 4 C, see also right hemisphere and inflated view in Supplemental Figure S7).

Spatial overlap with probabilistic cytoarchitectonic atlas labels
Vertex-wise comparison of spatial overlap on Colin27's cortical surface reconstruction between eleven ex vivo probabilistic Brodmann areas ( Fischl et al., 2007 ) and matching regions in the digitized Brodmann surface atlas showed good to very good spatial overlap (average 0.79, std 0.15), indicating that the majority of the vertices belonging to the ex vivo labels fall within the region labels of the digitized Brodmann atlas (Supplemental Table S8).
Furthermore, voxel-wise analysis of region label overlap between the digitized atlas collection and probabilistic region labels included in the Julich-Brain cytoarchitectural atlas ( Amunts et al., 2020 ) resulted for each atlas, for each region in a list of most overlapping Julich-Brain labels (Supplemental Table S9). Overall, the largest overlap was observed in primary visual cortex (with Visual_hOc1_l, mean overlap 90.83%, std 9.22), while e.g. fronto-opercular regions (where the Julich Brain atlas is notably more granular) showed a lower degree of overlap (e.g. overlap with Broca_44_l on average 52.42%, std 30.30, with Brodmann and Kleist at 77%) (see Supplemental Table S9).

Spatial overlap with NeuroSynth functional annotations
The digitized atlases were applied to a pipeline overlaying each atlas region with 942 functional annotations from the NeuroSynth database of voxel-wise involvement in brain function. In each atlas, for each brain region, this resulted in a ranked top 10 of NeuroSynth terms with the highest involvement scores, giving an indication of putative fMRI-based brain functions and/or disorders localized in that brain region. The top 10 terms show the greatest consistency across atlases in areas involved and low (yellow, e.g. superior frontal cortex) consistency. C illustrates the spatial consistency of region boundary placement across atlases, showing (similar to the region consistency in B) highest consistency at borders of primary cortical areas, and more variable boundary placement in higher order areas of the cortex. The right hemisphere showed highly similar vertex-wise consistency maps (see Supplemental Figure S6 and S7). For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) in primary or unimodal processing (e.g. the visual or motor system, areas with also highly consistent spatial mapping across atlases), while for example prefrontal cortex shows much more atlas-specific top 10 functional terms (Supplemental Table S10).

Discussion
We present an atlas collection containing digital reconstructions of cortical atlases of pioneering neuroanatomists Campbell, Smith, Brodmann, Von Economo, Flechsig and Kleist. The atlases describe regional boundaries based on microscale histological information, including regional cytoarchitecture, myeloarchitecture, myelogenesis, main cortical types, and/or loss-of-function lesion mappings. Digital reconstructions of classic brain atlases are a potentially important resource for neuroimaging studies, widening the range of modalities for which cortexwide atlases are available. The atlases were evaluated in terms of spatial consistency with each other, as well as with external datasets on structural and functional annotation of the cortex. Within the atlas collection, the overall across-atlas spatial consistency was examined, highlighting zones of high consistency as well as specificity in region definition. Evaluating region overlap with the cytoarchitectural Brodmann labels in surface space ( Fischl et al., 2007 ) and the benchmark probabilistic Julich-Brain cytoarchitectural atlas in MNI space ( Amunts et al., 2020 ) shows good spatial label overlap. Functional validation of the atlases with NeuroSynth ( Yarkoni et al., 2011 ) annotations illustrates the potential functional involvement of atlas regions.
A number of remarks regarding the digital reconstruction of the atlases can be made. First, brain atlases were transformed from 2D drawings into a 3D digital model following the original documentation of the respective authors. The 2D-to-3D atlas translation process is inherently accompanied by a number of challenges. Some regional boundaries were described in more detail than others, which could influence the detail in which the manual delineation of cortical areas onto the Colin27 template could be performed. For example, some region borders (e.g. region 6 in the Flechsig atlas ( Flechsig, 1920 )) are described to lie 'at' a certain cortical sulcus. Others are depicted in a figure showing a region situated on the outer surface of the brain and delimited by a sulcus, with no additional textual clarification (e.g. region 21 in the Flechsig atlas ( Flechsig, 1920 )). In practical terms in the digitization process, such descriptions translated to a border placed at the floor of that sulcus. In contrast, of other region borders (such as e.g. area 4 in the Brodmann atlas ( Brodmann, 1909 ) more information was available that could be used for the digitization process.
Second, the original atlases were in most cases constructed based on the examination of multiple brains, but we were not always able to find the exact detailed demographics of the subjects used in the original documents (information that we were able to find is presented in Table 1 ), limiting sample generalizability. Third, the original atlases were constructed as an averaged regional subdivision across multiple brains, resulting in static group-level cortical parcellations presented to the reader on a single cortical surface depiction. As such, as a limitation, the presented classic atlases do not provide information on individual variation, data that is available in probabilistic cortical mappings such as the Julich Brain atlas ( Amunts et al., 2020 ). A fourth limitation is that the classic atlases were digitized on the Colin27 template, a standard brain template generated based on 27 structural scans of a single individual. We evaluated the effect of atlas construction on an increasing number of individual segmentations on atlas segmentation consistency (see Supplemental Information, and Supplemental Figure S12), and though atlas consistency is increased with incorporation of more individual segmentations, single-segmentation atlases performed reasonably well. Still, the single-template basis of the presented atlas collection is a potential confounding factor that should be kept in mind when implementing these atlases in analysis pipelines. Fifth, we were not able to find information on lateralization in the original documentation, with no clear information on whether segmentations were performed based on the left, right or both hemispheres. As a result, the presented digital versions of the atlases also do not contain information on potential lateralization or left-right asymmetries of region boundaries. The digital reconstructions do include area segmentations for both the left and right hemisphere to accommodate neuroimaging analyses in which bihemispheric measurements are available, this was achieved by applying the documented information to both hemispheres. Sixth, the Neu-roSynth database ( Yarkoni et al., 2011 ) presents a unique resource of fMRI terms and associated spatial activation patterns. The meta-analysis data mining approach of the database can provide a general overview of the meta-activation maps of functional terms, and the overlap with the digital atlases as presented in this paper should thus be taken as such, providing a general indication of possible functional involvement of the regions indicated in the presented atlases. It further has been mentioned that MNI space volumetric across-subject registration methods tend to have relatively lower spatial precision than surface-based alignment methods ( Coalson et al., 2018 ). The voxel-wise evaluation analyses should be taken as a general indication of functional and cytoarchitectural co-localization with digitized atlas regions. Future studies employing more tailored, higher resolution surface-based maps (e.g. Glasser et al., 2016 ) would be of great interest and could provide more spatial specificity of effects.
Within the neuroimaging field, a wide range of cortical atlases are being used, providing means to parcellate the cortex based on all sorts of modalities and in varying resolutions ( Dickie et al., 2017 ;Evans et al., 2012 ). With an ongoing presentation of functional and structural atlases of the human cortex based on neuroimaging-derived features, classical atlases embody different, complementary aspects of micro-structural cortical organization and digital reconstructions of these atlases may therefore be of great added value to the MRI community. The variation in microstructural properties of cortical areas captured in these classical atlases -together with other properties such as neurotransmitter receptor densities, the level of interplay between cells in a cortical column and level of connectivity between brain areas-shape the cortical substrate underlying the function of an area, and the role it plays in the macroscale organization of the human cerebral cortex commonly evaluated using in vivo neuroimaging techniques ( Caspers et al., 2013 ). Integrating information on regional cytoarchitecture, myeloarchitecture and cognitive functionality in MR processing pipelines might therefore provide a more extensive context to research questions in neuroimaging research ( Dickie et al., 2017 ;Evans et al., 2012 ). Incorporation of the presented collection of digitized cortical atlases into neuroimaging workflows may further aid in consolidating neuroimaging findings across atlases ( Madan and Kensinger, 2018 ;Shafiei et al., 2020 ), facilitate anatomy-guided node selection for connectome analyses, provide ROIs that can be used to test anatomically specified hypotheses ( Eickhoff et al., 2006 ) and serve as new parcellation schemes for neuroimaging studies that focus on bridging the gap between micro-and macroscale neuroscience (

Data availability statement
All digital brain atlases are publicly available and free to use under a Creative Commons Attribution-NonCommercial-Share A like 4.0 International License.
Atlas packages are available for download at www.connectomelab.nl/ (downloads). The Colin27 and ICBM152 brain templates can be found at www.bic.mni.mcgill.ca . The data of the Human Connectome Project data is freely available at www.humanconnectome.org . The funding agencies did not have any influence on data acquisition, analysis and/or report of the data. All authors report no biomedical financial interests or potential conflicts of interest.