2D histology meets 3D topology: Cytoarchitectonic brain mapping with Graph Neural Networks

Cytoarchitecture describes the spatial organization of neuronal cells in the brain, including their arrangement into layers and columns with respect to cell density, orientation, or presence of certain cell types. It allows to segregate the brain into cortical areas and subcortical nuclei, links structure with connectivity and function, and provides a microstructural reference for human brain atlases. Mapping boundaries between areas requires to scan histological sections at microscopic resolution. While recent high-throughput scanners allow to scan a complete human brain in the order of a year, it is practically impossible to delineate regions at the same pace using the established gold standard method. Researchers have recently addressed cytoarchitectonic mapping of cortical regions with deep neural networks, relying on image patches from individual 2D sections for classification. However, the 3D context, which is needed to disambiguate complex or obliquely cut brain regions, is not taken into account. In this work, we combine 2D histology with 3D topology by reformulating the mapping task as a node classification problem on an approximate 3D midsurface mesh through the isocortex. We extract deep features from cortical patches in 2D histological sections which are descriptive of cytoarchitecture, and assign them to the corresponding nodes on the 3D mesh to construct a large attributed graph. By solving the brain mapping problem on this graph using graph neural networks, we obtain significantly improved classification results. The proposed framework lends itself nicely to integration of additional neuroanatomical priors for mapping.


Introduction
Cytoarchitectonic areas are characterized by a distinct spatial organization of neuronal cells in the brain, including their arrangement into layers and columns, as well as density, orientation, or presence of certain types of cells. As indicators for connectivity and functional modules, they provide an important microstructural reference for human brain atlases [2]. Mapping the borders of cytoarchitectonic areas requires to analyze brain sections at microscopic resolution. To capture the human brain's considerable variability, it needs to be performed in many different brain samples. [2] created a comprehensive probabilistic cytoarchitectonic atlas based on delineations in serial histological sections stained for cell bodies in a sample of ten brains, using image analysis and statistical tools to identify architectonic borders [22]. These resulting maps are aggregated in a 3D reference space at the millimeter scale. Although today's high-throughput scanners allow to digitize complete human brains at micrometer resolution, it is practically impossible to scale this approach for doing delineations in all sections of a human brain, which may have, in dependence of the size, approx. 6000-8000 sections. This motivates the development of automatic mapping algorithms.
Recent work formulated cytoarchitectonic mapping of cortical regions as an image segmentation [23,24] or classification [21] problem, which can be approached with deep convolutional networks. These methods process each 2D section individually, responding to the lack of routine workflows for computing a precise 3D reconstruction from histology en par with high throughput imaging. In fact, as of today only one 3D reconstruction of a whole human brain from histology is available [1], with a spatial resolution of 20 micrometer isotropic. The 3D topology of a brain is thus not used by these models, although it is key to disambiguate complex or obliquely cut brain regions. Other authors suggested to incorporate inter-slice information using 3D convolutions [19] or recurrent networks [3], which again requires a precise 3D reconstruction.
Here we present a novel paradigm for cytoarchitectonic brain mapping which overcomes the above restrictions. The basic idea is to reformulate the mapping task as a node classification problem on the approximate 3D midsurface mesh through the isocortex. Since the mesh is not assumed to be precise, it can be derived from a simple linear reconstruction which is straightforward to compute from a histological image stack, and can handle a significant amount of missing sections. Building on [21], we extract deep features from 2D cortical patches at microscopic resolution using convolutional neural networks (CNNs) that were trained with a contrastive learning task to encode cytoarchitectonic characteristics. These features are then assigned to the corresponding nodes on the reconstructed surface mesh to construct a large, attributed graph. The brain mapping problem is solved on this graph using graph neural networks (GNNs) [11,8].
We make the following contributions: 1) We introduce a novel deep learning approach for cytoarchitecture classification in large stacks of whole brain sections which integrates highresolution 2D texture features with global 3D topology. 2) The approach outperforms the state of the art on DL-based cytoarchitectonic mapping on a dataset of histological sections from eight postmortem human brains. 3) The approach allows flexible integration of neuroanatomical priors into the mapping problem, which further boosts classification performance.

Methods
The proposed framework consists of three components: 1) A linear 3D reconstruction of the histological stack to compute an approximate midsurface mesh of the isocortex, which is interpreted as a graph. 2) A CNN trained with contrastive learning which extracts cortical features from the histological sections, serving as node attributes in the graph. 3) A GNN to label each graph node with a cytoarchitectonic area, exploiting both high-resolution texture and approximate 3D topology encoded in the graph.
Construction of the approximate midsurface graph A coarse segmentation of each section into background, white matter and cortical gray matter is performed using morphological active contours [17]. Histological brain sections (Fig. 1, left) are 3D reconstructed by rigid alignment of adjacent sections using feature based registration as described in [6]. From the segmented linear 3D reconstruction, a 3D cortical midsurface mesh is extracted by solving Laplace's equation inside the cortical mantle [13] using the highres-cortex package 1 and extracting the 0.5 isosurface from the resulting volumetric Laplacian fields using marching cubes [14]. In a manual postprocessing stage, hemispheres are split, and small disconnected components as well as cerebellum and brainstem are removed from the mesh. Finally, isotropic explicit remeshing is applied using meshlab [4] to unify edge lengths in the mesh to ≈ 300µm, thereby reducing the number of triangles and ensuring that connections between vertices represent comparable distances. Each resulting mesh (Fig. 1, right) is then interpreted as a graph G = (V, E) with nodes V and edges E. Computing node level texture features Each node u ∈ V in the resulting graph can be uniquely identified with an image coordinate p u in a 2D histological section by inverting the rigid transformation applied for 3D reconstruction. We apply an encoder CNN f to extract a cytoarchitectonic feature embedding h u from an image patch x centered at each coordinate p u . The encoder f is trained using the contrastive approach proposed in [21] to produce an embedding that maps image patches from the same brain area to similar feature vectors, and image patches from different brain areas to dissimilar feature vectors. Given an image patch x, f produces a lower dimensional vector h = f (x) ∈ R De , which is passed through a projection network g(h) =z ∈ R Dp and normalized as z i =z i / z i 2 . Given a minibatch of N image patches x i with corresponding labels y i (i = 1, . . . , N ), the contrastive loss is optimized during training, where I is the indicator function, τ ∈ R is a temperature scaling parameter and N y i is the number of batch elements with the same label as x i . The projection network g is discarded after training.
Graph neural networks We consider two GNN architectures to integrate spatial relationships of image patches into the classification model: GraphSAGE [8] and Graph Attention Network (GAT) [25]. Both architectures are suitable for inductive tasks, so that trained networks can be applied to previously unseen nodes or graphs. Due to common memory constraints, training on full graphs is not possible given the level of detail of our models. We therefore adapt the neighborhood sampling scheme proposed in [8]: Given a node u ∈ V and a GNN with K layers, we sample its K-hop neighborhood N K (u) and use the resulting subgraph as input for the GNN. This way, information from the immediate neighborhood of u can be propagated through the GNN to classify u.
Integrating anatomical priors Previous work [23] showed that the integration of prior anatomical knowledge -there given in the form of probabilistic cytoarchitectonic maps [2] -can improve classification performance. Inspired by this finding, we adapt the registration workflow from [2] to project probabilistic maps of 152 brain areas from the MNI Colin27 reference space onto the histological sections, and assign an additional vector h P u to each node u. Each dimension of h P u encodes the probability to belong to a particular region of the Julich-Brain atlas [2]. We further consider spatial locations in Colin27 space as another anatomical prior by assigning a 3D coordinate vector h C u to each node u, again using the workflow from [2].

Experiments and Results
We systematically evaluate all components of the proposed approach: 1) We compare different encoder architectures for contrastive feature learning from images patches. 2) Using the best performing encoder architecture, we compare the performance of GNNs with SAGE and GAT architectures. 3) We investigate the additional benefit of adding neuroanatomical priors for classification.
Dataset We use an in-house dataset containing images of 1860 histological sections from 7 human postmortem brains acquired from the body donor programs of the Anatomical Institute of the University of Düsseldorf (Germany), and corresponding annotations of 113 cytoarchitectonic cortical areas. Sections have an approximate thickness of 20µm, were marked for neuronal cell bodies using a modified Merker stain [18], and imaged at a resolution of 1µm with a lightmicroscopic scanner (TissueScope, Huron Digital Pathology Inc.). 80% of the sections are used for training (X tr ), remaining sections for testing (X te ). Transferability to unseen brains is evaluated on 325 sections from an 8th brain (X un ). Graph nodes V are split into V tr , V te and V un , containing nodes with corresponding points p u in X tr , X te and X un , respectively. Performance is evaluated using macro-F1 score on nodes V te and V un . The study requires no ethical approvals.

Performance of different feature encoders
The encoder network f is trained on 1200 image patches per class sampled from X tr , oversampling small areas for class balancing. Data augmentation steps include random pixel intensity augmentation, rotation, mirroring, translation, blurring and sharpening with parameters detailed in [21]. Training is performed for 150 epochs using LARS optimizer [26] with constant learning rate 0.01 * N/128, batch size N = 4096 (N = 2048 for res[w,l] due to memory constraints) and τ = 0.07. We consider the two network architectures listed in  [10], ReLU) to classify brain areas from features extracted by each considered model f for all nodes in V tr . MLP is trained for 100 epochs using SGD with Nesterov momentum 0.9, constant learning rate 0.001 * N/4096, batch size N = 16384 and cross-entropy loss. We apply dropout with probability 0.5 and 0.25 to the input and hidden layers, respectively. Table 1: Network architecture of base and res models. res uses pre-activation residual connections as in ResNet18 [9]. base uses no padding, res uses zero-padding (CONV (k  We found that res models perform better than base models (Table 2 left). Wider networks and a larger input size improve performance on V te . On V un , increased input size (res[l]) outperforms wider networks (res[w]). Best results are obtained by combining both (res[w,l]).
Graph neural networks We create one graph per brain hemisphere, resulting in 16 graphs overall (1.2 million nodes with node degree 6, on average). Training nodes are sampled from V tr , inversely proportional to their class frequency to account for class imbalance. Node features are computed using res[w,l]. Node labels are assigned using annotations of cytoarchitectonic areas. We consider two basic GNN architectures: SAGE [3] consists of three GraphSAGE layers with mean aggregation [8], while GAT [3] consists of three GAT layers [25]. Following [15,16] [8] with a neighborhood size of three. For GAT models, we sample full neighborhoods, but apply dropout with probability 0.5 to attention coefficients as proposed in [25]. Both sampling methods aim to make models robust against missing nodes or edges.
All GNN architectures obtain significantly higher scores than MLP on both V te and V un ( Table 2, center). SAGE[5+r] achieves overall best scores, with an increase by 19.61 points and 1.23 points on V te and V un compared to MLP, respectively. Using the same number of layers and residual connections, SAGE outperforms GAT in almost all cases. Incorporating prior knowledge We examine classification performance under all possible combinations of using cytoarchitectonic features extracted from images (CY), weights from probabilistic cytoarchitectonic maps (PM) [2], and canonical 3D locations (CO) at each node. Dropout with probability 0.5 and additive Gaussian white noise with standard deviation 0.05 is applied to h P u and h C u , respectively. h P u and h C u are each projected by a fully connected layer (256 units, BN, ReLU) before concatenation.
Incorporating both PM and CO with the texture features CY improves scores by 9.75 points on V te and 12.45 points on V un (Table 2, right; Fig. 2). Using them in isolation, CO leads to slightly better performance than PM. Scores on V te drop significantly when removing texture features CY. Best scores on V un are obtained using only canonical coordinates, without texture features.

Discussion & conclusion
We presented a graph neural network approach for cytoarchitectonic classification of cortical image patches in microscopic scans of human brain sections, which combines 2D texture features with topological information from an approximate 3D surface reconstruction. While [5] employed GNNs for surface parcellation of Broca's areas 44 and 45 using MRI-based features, the present work is the first to our knowledge that integrates deep texture features from histology with 3D topology in a graph framework to label a large number of highly different brain areas across several brains. The proposed method significantly outperforms previous methods that operate only on individual 2D sections: The best trained model increases classification scores by 161%/193% on V te /V un wrt. the recent baseline model base representing the work from [21]. Our experiments further suggest that deeper and wider architectures for the CNN encoder f benefit performance, motivating more systematic architecture search in the future. Larger input sizes in the encoders res[l] and res[w,l] showed additional positive impact by enabling models to capture the entire width of the isocortex (2mm-4mm). SAGE outperformed GAT, suggesting that mean aggregation (SAGE) might be better suited than self-attention (GAT) for our use case. The presented framework allows straightforward integration of anatomical priors, and the results indicate that this might be a crucial strategy for optimizing cytoarchitectonic mapping with deep networks. Unfortunately, transferability of learned features to completely unseen brains still seems to be limited, as indicated by the performance gaps between V te and V un and thus confirming observations in [21]. In the future, we plan to perform more rigorous architecture search for the encoders, study more anatomical priors, and investigate the reasons underlying the reduced performance on unseen brain samples.