Histological data were obtained from the Zagreb Neuroembryological Collection24. Samples used in this study were taken from the prefrontal cortex of two brains (brain 1 55 years old female, post-mortem delay 24h; brain 2 age not available; male, post-mortem delay 4h). Sections were taken from the dorsal and ventral part of the typical six-layered homotypic isocortex of the prefrontal cortex3,4. Brains were fixed in \(4\%\) PFA for two weeks. Following sampling, sections were dehydrated through a series of ethanol and embedded in paraffin. Sections were cut using rotating microtome at thicknesses of \(10\;\upmu \hbox {m}\) and \(20\;\upmu \hbox {m}\). The tissue was stained using the NeuN immunohistochemistry method according to standard protocol25. NeuN is an RNA-binding nuclear protein, derived from the RBFOX3 gene, which regulates alternative splicing in neurons and is expressed explicitly in all neurons of used tissue specimens. In the experiments, \(10\;\upmu \hbox {m}\) and \(20\;\upmu \hbox {m}\) sections were used in order to test whether tissue thickness would impact the results. Histological sections were digitized using the Hamamatsu Nanozoomer 2.0 scanner (Hamamatsu Photonics, Japan) at 40x magnification, corresponding to \(0.226\mu \hbox {m/pixel}\) resolution. Example histological sections in Fig. 1 show varying neuronal morphology and cellular distribution across the cortical layers. Computational experiments were performed using custom scripts written in Python 3.8 and standard publicly available libraries.
Histological slices stained with NeuN immunohistochemistry method showing varying neuronal morphology and cellular distribution across the layers of the cortex. The left image is taken from the dorsolateral part of the prefrontal cortex (thickness \(10\;\upmu \hbox {m}\)) and the right image was taken from the orbital part of the prefrontal cortex (thickness \(20\;\upmu \hbox {m}\)). Scale bar \(100\;\upmu \hbox {m}\).
Neuron-level features
In manual delineation, the density and size of neurons are the most important characterizations of the laminar structure. Based on anatomical descriptions and kernel density estimates, three populations of similar densities are assumed (layers II and IV as dense, layers III, V, and VI as average, and Layer I and white matter as sparse), two populations of similar sizes (layers III, V and VI containing on average larger neurons, and layers I, II, IV and white matter containing on average smaller neurons). By plotting a histogram of neuron densities and sizes, one can observe that the features express a multimodal distribution, which can be separated using the minimization of intraclass variance26. Figure 2 shows a visualization of separating the neuron populations across the histological section, revealing the laminar structure. In contrast to the classical pixel-based approach to cortical layer segmentation, we used neuron-level tissue descriptors to characterize and examine the underlying tissue properties. We develop several feature classes that describe each neuron in the tissue and use a machine learning model to determine the layer of individual neurons. By classifying all neurons in the tissue, we obtain the parcellation of the laminar structure. To the best knowledge of the authors, this is the first bottom-up approach in the analysis of brain cytoarchitectonics that builds from the cellular level and infers about larger structures based on morphological and textural features of individual neurons. Here, we discuss the development of those features, as well as some rationale behind the choices made in their development and selection.
The first step in obtaining the neuron characterization is the segmentation of neurons from the background tissue. The segmentations were obtained using automated methods27,28 which use grayscale-guided watershed on anisotropically diffused images to separate neurons, rather than often used distance maps obtained from the grey-level threshold, providing a binary image of segmented, non-overlapping neuron areas. As the goal here is to create consistent results across the tissue, other segmentation methods may be used as well, such as a recently proposed instance segmentation via contour proposals29, especially for different staining methods. This step yielded the locations and segmentations of neurons, from which other neuronal characteristics are developed.
Secondly, neuron segmentations were analyzed using ImageJ particle analysis pipeline30. An overlap of binary segmentations and the original image was made, and ImageJ’s function analyze particles was used to produce measurements of neurons’ bodies. Those were the area, perimeter, circularity, roundness, and Feret’s diameter as well as the mean, median, skewness, and kurtosis of the gray values. More details on particle measurement can be found in ImageJ’s documentation31. These features form the basis for investigations in brain microanatomy, as they are often, although not at this level of precision, perceived by the eye of neuroanatomists. By this, we obtain the first neuronal characteristics, which may be visualized to reveal patterns of their appearance across the cortical layers, as shown in Fig. 2. It should be noted that values based on image intensity were not used in the further analysis as it was concluded that these were not usable generally, as they may be heavily influenced by uneven staining across the section and exhibit different values for different staining procedures. These simple measurements do not possess a discriminative power to create clear classifications of neurons within the layers. Therefore, richer descriptors that incorporate neuron neighborhoods were computed, as described below.
By visualizing the basic neuronal characteristic, the area, and local density, the laminar structure is revealed, as these features exhibit multimodal distributions. Left: Three types of neurons were distinguished by cell density in their surroundings: very sparse (blue), sparse (red), and dense (green). Middle: Larger (red) and smaller (green) neurons. Right: Average size of neighboring neurons is a feature obtained from the previous two: nearest neighboring neurons and their areas, which can also be used to facilitate layer prediction. Here, thresholds that separate the distributions were obtained using minimization of intraclass variance26.
It is worth mentioning that density-based clustering algorithms are often used to segment the areas of similar point densities32,33,34, which may be brought into relation with cortical layers having a roughly uniform density within each layer. However, it seems that the neuron distribution in the cortex is such that their intrinsic structure may not be clustered by a single set of global density parameters, as used in clustering methods. Nevertheless, these methods provided insight into some cortical properties. Meaningful clusters were created when considering neurons within the radius between \(100\;\upmu \hbox {m}\) and \(300\;\upmu \hbox {m}\), containing between 300 and 800 neurons. This lead to the conclusion that the changing nature of neuron distribution in the brain is best characterized when performing measurements in this range. This range approximately also corresponds to the biological limits of interlayer distances. It is important to emphasize that although choosing a predefined radius or a number of neighbors may seem equivalent, analysis of nearest neighbors is preferred over the fixed radius approach. A predefined radius may be interpreted differently, depending on the image resolution. The specified range around a neuron allows for a more detailed and precise analysis of the microstructure of the tissue in that specific area by balancing between obtaining enough information to capture local tissue properties and not being confounded by reaching too far from the neuron toward other layers and incorporating information which is not in the neuron’s vicinity and is, therefore, less relevant for neuron’s phenotype. Also, if a fixed maximum number of neighbors for each neuron is used, efficient data structures like kd-trees35,36 may be precomputed. Considering the large number of neurons found in a histological section, efficiency may be of critical importance.
To measure properties of neurons’ neighborhoods, nearest k neighbors were considered, for \(k \in [50,100,250,500,1000]\). The distances to a neuron’s k-th nearest neighbor were used as a feature, as well as their mean, max, min, skewness, kurtosis, and entropy. Basic measures of individual neurons were computed in a similar fashion to produce, for instance, the average area of neighboring 100 neurons, as shown at the right in Fig. 2. A convex hull of neurons’ k-neighbours gives information about the area around a neuron and a number of its neighbors and is described using hull area, perimeter, average nearest distance for neurons found in the hull, and standard deviation of nearest distances. Dispersion of neurons may be quantified using nearest neighbor index (NNI), a measure that describes whether points follow usually subjective patterns of regular, clustered, or random distribution. The NNI measures the distance between each point and its nearest neighbor’s location. All the nearest neighbor distances are averaged, and if the average distance is less than the average for a random distribution, the distribution of the features being analyzed is considered clustered. If the average distance is greater than a random distribution, the features are considered regularly dispersed. The index is expressed as the ratio of the mean observed distance divided by the expected distance, which is based on a random distribution with the same number of points covering the same total area,
$$\begin{aligned} NNI_i = \displaystyle \frac{\frac{1}{n} \sum _{j=1}^n d(i,j)}{0.5 \sqrt{HullArea(i)/n}}. \end{aligned}$$
(1)
Neurons in all layers except layer I and white matter tend more towards uniformly dispersed distribution, especially neurons of layer IV which tend more towards random distribution.
Depending on its position in the cortex, a neuron may be placed more toward the middle or more toward the edge of its layer. The computation of properties of its neighborhood may be confounded by reaching into adjacent layers and using neurons with different properties for the computation of statistics. To identify this case, measurements may be taken only from neurons found within the range of angle, or slices. Features measured in several directions can identify border neurons and changes in neuronal properties in different directions. Slices may be regarded as measurement units reaching from a single neuron, each unit representing a population of neighboring neurons found in a given direction from the central neuron. The relationship of different populations within an area has been extensively studied in the frame of biological diversity of species, landscapes and other37,38. Considering the neurons in a slice as members of a single species, and the k neighbors of a neuron as the population of all species in their habitat, biodiversity measures evaluate the relationship between the species. In this context, the number of slices is the number of different species, or richness, and the relative abundance of the different species in an area as evenness. The two most often used such measures are the Shannon index39 and Simpson index40. The Shannon index gives a quantitative measure of the uncertainty in predicting the species of an individual chosen randomly from the population. The Simpson index measures the probability that the two individuals who are randomly chosen (with replacement) from the total population will be of the same species.
$$\begin{aligned} Shannon = – \sum _{i=1}^{R} p_i \ln p_i = \ln \left( \frac{1}{\prod _{i=1}^{R} p_i^{p_i}} \right) , \quad Simpson = \sum _{i=1}^{R} p_i^2, \end{aligned}$$
(2)
where R is the number of different species or, here, slices, and \(p_i\) is the proportion of species of the ith type in the population or proportion of neurons in ith slice to the number of the neurons in k-neighbourhood. If all slices have an equal number of neurons, \(p_i\) values equal 1/R, and the Shannon index takes the maximum value of \(\ln R\). If the numbers are unequal, the weighted geometric mean of the \(p_i\) values is larger, which results in the index having smaller values. The index equals zero if the neurons from only one slice are present since there is no uncertainty in predicting the slice they are in. The index gives information about the relation between the number of types and the presence of the dominant type. The mean proportional abundance of the slices increases with decreasing number of slices and with the increasing abundance of the slice with the largest number of neurons, the index obtains small values in regions of high diversity like neurons on borders between the layers, thin layers, and especially layer I neurons. The index is large in homogeneous areas like the middle of layer III, where slices reaching from a neuron remain in the area of the layer.
Experimental subjects statement
All specimens were collected during regular autopsies at pathology departments of the University of Zagreb, School of Medicine, approved by the Ethics Committee of the University of Zagreb, School of Medicine and in accordance with the Declaration of Helsinki, and informed consent was obtained from the next of kin.
