BasNet: Внимание, У-Net-автоматизированная сегментация аксонов в Bielschowsky Silver-Stained Histology.

Текст сохранён. Часть иллюстраций пока доступна только в первоисточнике; загрузка будет проверена повторно.

Introduction

Precise, reproducible quantification of axonal density is essential across neuroscience and neuropathology applications because axons constitute a major element of structural integrity and functional capacity of neuronal tissue (Kevenaar & Hoogenraad, 2015 ; Stadelmann et al., 2019 ). Quantitative measures of axonal loss or remodeling are central to the study of disease mechanisms combining histopathology and imaging (Gkotsoulias et al., 2025 ) and to the objective assessment of progression and therapeutic effects in conditions such as multiple sclerosis (MS) and other neurodegenerative disorders.

Silver impregnation methods - most commonly the Bielschowsky silver stain - remain widely used to visualize axons and neurofibrils in clinical and research histology (Bielschowsky, 1902 ; Nestor, 2008 ). However, converting silver-stained images into reliable quantitative maps remains challenging. Manual counting and stereological sampling are labor-intensive and scale poorly, while surrogate measures such as optical density and line-intersection counts are confounded by staining and sampling choices: optical density depends on staining quality and also captures non-axonal Bielschowsky-stained tissue, while line-intersection counts vary with arbitrary line placement and local fiber orientation across large tiles or whole-slide images (Fathy et al., 2022 ). In practice, these limitations reduce throughput and can introduce systematic bias that complicates comparisons across cohorts and studies.

Recent advances in convolutional neural networks have produced robust tools for axon and myelin segmentation in a variety of imaging modalities - for example AxonDeepSeg (Zaimi et al., 2018 ) and DeepACSON (Abdollahzadeh et al., 2021 ) - demonstrating that deep learning can deliver accurate, high-throughput morphometric analyses of neural tissue. Nevertheless, most existing deep learning pipelines were developed and validated on modalities such as electron microscopy, immunofluorescence, or myelin-specific stains, and to our knowledge no models optimized for brightfield silver-stained histology such as Bielschowsky exist. This modality is widely used in research and the clinical context, but presents distinct challenges (heterogeneous contrast, strong background staining) that can degrade the performance of models trained on other modalities.

Here we present an attention U-Net-based segmentation model (Ronneberger et al., 2015 ; Oktay et al., 2018 ) trained on manual annotations of axons and neurites sampled from multiple brain regions and from demyelinated lesions. The model produces dense, tile-level and whole-slide axonal segmentation maps that enable automated, reproducible axon density estimation and permit downstream morphometric analyses (e.g., fiber density maps and fiber-orientation/dispersion metrics). Compared with manual stereology or intensity-based surrogates, our approach increases throughput, removes arbitrary sampling decisions, and yields spatially resolved, directly interpretable axon maps that facilitate within-slide and between-cohort comparisons. We also provide a self-contained software package for straightforward deployment on tiled images and whole-slide scans, enabling integration into routine pathology-informatics workflows.

Contributions

A robust attention U-Net trained on manually segmented Bielschowsky silver-stained histology, validated across multiple tissue classes (white matter, gray matter, demyelinated lesions, and background).

Demonstration of downstream utility: automated axon density maps and fiber-orientation metrics derived from segmentations.

A ready-to-deploy software application for tile-level and whole-slide processing to enable reproducible, high-throughput axon quantification.

Methods

Tissue Preparation and Histological Processing

The tissue samples used in this study were dissected from four brains from patients with MS, provided by the German MS Brain Bank of the Competence Network Multiple Sclerosis. This study was approved by the ethical review committee of the University Medical Center Göttingen. Representative samples of gray matter, white matter, and MS-lesion areas were manually chosen in each tissue sample. From within these, tiles sized 1024² pixels (pixel size 0.273 μm) were selected, and contrast enhancement was performed using Contrast Limited Adaptive Histogram Equalization (CLAHE; Zuiderveld, 1994 ) (clip limit: 4.0, grid size: 8 pixels). Manual ground truth segmentation of axons and dendrites on the resulting images was performed by two medical doctors. The two raters divided the annotation workload, so that each tile was segmented once by a single rater rather than by both (see Sect. 2.4 ). Cell bodies of all cell types were omitted from the segmentation. The same process was used for the validation and test sets (512² pixels). The reduced tile size was chosen as a practical tradeoff, balancing segmentation workload and coverage of tissue variability across the dataset. We deliberately selected tiles from diverse anatomical regions, tissue types, and WSIs to maximize biological and technical heterogeneity. A standardized 256 × 256 patch extraction ensures that smaller validation/test tiles do not introduce distributional shift and that the model observes identical spatial contexts regardless of source tile size. The finished ground truth masks were binary PNG images with values 0 for background and 255 for axons.

Network Architecture

We developed an automated segmentation pipeline using an Attention U-Net architecture for binary classification of axonal structures in Bielschowsky-stained tissue sections. The network architecture followes the encoder-decoder structure originally proposed by Ronneberger et al. ( 2015 ) and was enhanced with attention gates as described by Oktay et al. ( 2018 ) to selectively emphasize task-relevant features while suppressing background information.

The encoder consisted of four hierarchical levels with progressive channel expansion (32, 64, 128, 256 channels) and spatial downsampling via max pooling operations. Each encoder level comprised two sequential convolutions, each followed by batch normalization and rectified linear unit (ReLU) activation. Dropout regularization ( p=0.3 ) was applied at encoder levels 2–4 to prevent overfitting. The bottleneck layer expanded to 512 channels at an spatial resolution.

The decoder mirrored the encoder structure with four upsampling levels, utilizing transposed convolutions with stride 2. At each decoder level, attention gates modulated the corresponding encoder features before concatenation with the upsampled decoder features. The attention mechanism computed spatial attention coefficients using the formulation:

where g represents the gating signal from the decoder, x denotes the skip-connection features from the encoder, and are convolutions for channel reduction, is a convolution to a single channel, and is the sigmoid activation function. The attention-weighted features were then computed as the element-wise product .

The final output layer consisted of a convolution producing single-channel logits for binary segmentation. The complete network contained approximately 7.85 million trainable parameters with input dimensions of pixels (RGB) and output dimensions of pixel (binary mask).

Loss Function and Optimization

To address the severe class imbalance inherent in axon segmentation (approximately 30% foreground vs. 70% background pixels), we employed Focal Tversky Loss (Abraham & Khan, 2018 ), which combines asymmetric penalty weighting with hard-example mining. The Focal Tversky Loss is defined as:

where TI represents the Tversky Index (Tversky, 1977 ):

with TP, FP, and FN denoting true positives, false positives, and false negatives, respectively. We set (false positive weight), (false negative weight), (focusing parameter), and (smoothing constant). The asymmetric weighting ( ) prioritized recall over precision to minimize missed axon pixels, while the quadratic focusing parameter ( ) up-weighted difficult examples, directing learning toward challenging boundary regions. The choice of β = 0.6 is empirically supported by a post-hoc parameter sweep (Sect. 3.6 ).

The network was optimized using AdamW (Loshchilov & Hutter, 2017 ) with the hyperparameter settings used in the released training configuration. Learning rate scheduling employed ReduceLROnPlateau with a reduction factor of 0.5 and patience of 5 epochs. Gradient clipping with a maximum norm of 1.0 was applied to stabilize training. Training continued for a maximum of 100 epochs with early stopping (patience: 5 epochs) based on validation loss.

Dataset Composition and Annotation Strategy

Manual annotation of axonal structures in Bielschowsky-stained tissue is extremely time-intensive, necessitating a pragmatic approach to dataset assembly. Two trained medical doctors annotated 33 tissue tiles from 26 distinct whole-slide images (WSIs) across four human brains (brains 1, 2, 3, and 4), capturing white matter (WM), gray matter (GM), background (BG), and demyelinated lesions (L) across multiple anatomical regions.

To make this labor-intensive task feasible, the two raters divided the annotation work rather than each annotating every tile: for model training and evaluation, every tile carried a single annotation from one rater, so the model never received two conflicting labels for the same image and no label fusion was performed. Before annotating, both raters agreed on a common segmentation protocol by jointly performing segmentations on several shared example tiles. The segmentations was mainly based on two rules: (i) segmenting only thin, dark silver-impregnated axonal fibers, excluding cell somata, and, (ii) where an axon passed through a soma, tracing it continuously through the soma rather than interrupting it. The first rater annotated all validation and test tiles and 13 of the 17 training tiles, while the second rater annotated the remaining 4 training tiles. The model was therefore trained and evaluated against a single, consistent annotation standard (all test tiles originate from the first rater) rather than a fusion of two rater styles.

WSI identifiers follow the format brain_block (e.g., 1_1, 1_2 represent different tissue sections from brain 1). The dataset was strategically partitioned to balance regional representation, maximize test set independence, and ensure consistent tile dimensions within validation and test sets:

Training set : 17 tiles from 12 WSIs across 3 brains (1, 2, 3).

Validation set : 8 tiles from 8 WSIs across 2 brains (3, 4).

Test set : 8 tiles from 8 WSIs across 4 brains (1, 2, 3, 4).

The dataset composition is summarized in Table 1 . Critically, 7 of 8 test WSIs (87.5%) were completely independent, never appearing in training or validation sets. Only one test WSI (2_3) was shared with the training set, but tiles represented spatially distinct tissue regions (WM in test set vs. GM in training). Similarly, only one WSI (3_13) contributed tiles to both training and validation sets, again with spatially non-overlapping regions (BG in validation set vs. GM in training). No WSIs were shared across all three data splits.

Regional composition was balanced across all sets to ensure model exposure to diverse tissue types during training and robust evaluation across anatomical contexts. Training tiles comprised 5 WM, 5 GM, 4 Lesion, and 3 BG regions. Validation and test sets maintained similar regional distributions (WM: 3, GM: 2–3, Lesion: 1–2, BG: 1 per set), enabling assessment of generalization across tissue categories.

Detailed WSI Assignments

Training WSIs (12) : 1_1 (3 tiles), 1_4 † (2 tiles), 2_1 (2 tiles), 2_2 (1 tile), 2_3 ‡ (1 tile), 3_1 (1 tile), 3_5 † (2 tiles), 3_6 (1 tile), 3_8 (1 tile), 3_12 (1 tile), 3_13 † (1 tile), 3_15 (1 tile).

Validation WSIs (8) : 3_2 * , 3_3 * , 3_9 * , 3_10 * , 3_11 * , 3_13 † , 3_16 * , 4_3 * (all 1 tile each).

Test WSIs (8) : 1_2 * , 1_3 * , 2_3 ‡ , 3_4 * , 3_7 * , 3_14 * , 4_1 * , 4_2 * (all 1 tile each).

* WSI unique to this split. † WSI shared between training and validation (3_13) or within training only (1_4, 3_5 multi-tile). ‡ WSI shared between training and test. Shared WSI always represent spatially distinct regions.

Patch Extraction Strategy

All tiles, regardless of original dimensions, were subdivided into standardized pixel patches, ensuring consistent spatial context for the model across all data splits. Patch extraction employed two complementary strategies: Systematic grid extraction with a 256-pixel stride ensuring complete tile coverage and uniform random sampling introducing positional variability.

The final patch-level dataset comprised:

Training : 17 tiles 86,330 patches after 5-fold augmentation.

Validation : 8 tiles 2,032 patches (no augmentation).

Test : 8 tiles 2,032 patches (no augmentation).

Data Leakage Prevention

To prevent data leakage through spatial autocorrelation, all patches extracted from a given source tile were assigned exclusively to a single data split. Adjacent tissue regions within the same tile exhibit structural similarity due to local cytoarchitecture, myeloarchitecture and staining characteristics, making random patch-level redistribution across splits methodologically inappropriate. Where WSIs contributed tiles to multiple data splits, tiles were spatially non-overlapping and sampled from distinct anatomical locations.

The high degree of test set independence (7/8 WSIs completely independent, including 2 WSIs from brain 4 never observed during training) enables robust assessment of model generalization. These design constraints rendered traditional k -fold cross-validation infeasible. Instead, we employed leave-one-WSI-out (LOWO) analysis on the training set to assess model robustness to exclusion of individual training WSIs.

Data Augmentation

Data augmentation was applied exclusively to the training set to simulate biological and technical variability, with five versions generated per original patch (1 original plus 4 augmented versions). Augmentation operations were applied using the Albumentations library (Buslaev et al., 2020 ). Spatial transformations were applied synchronously to images and masks, whereas intensity transformations were applied to images only.

Performance Evaluation

Model performance was assessed using multiple metrics computed on binarized predictions (threshold: 0.5). Primary metrics included the Dice coefficient (F1-score) and Intersection over Union (IoU, Jaccard Index). Secondary metrics included precision and recall. The Dice coefficient was calculated as:

and IoU as:

where X and Y represent predicted and ground truth masks, respectively, and prevents division by zero.

Interrater Agreement Evaluation

To directly quantify annotation subjectivity, two independent medical-doctor raters each segmented the same two image tiles. The resulting masks were compared pairwise with each other and with the model prediction using Dice coefficient and IoU. These two tiles — the test tile 1_3 and the validation tile 4_3 — are the only tiles annotated independently by both raters. For model training and evaluation only the first rater’s annotation of each was used, consistent with the rest of the validation and test set, which the first rater also labelled. The second rater’s annotations were made solely for this inter-rater comparison and were not used during the model’s training and validation.

Leave-One-WSI-Out Analysis

In each of 12 folds, one unique WSI was held out from training and the remaining WSIs were used to train a separate model. Validation and test set were held constant. Performance was then evaluated on the fixed validation set using Dice coefficient, IoU, precision, and recall. Each LOWO model was trained from scratch using the same architecture, loss function, optimizer, and early stopping configuration as the full model. Online Resource 1 summarizes the training loss curves for all folds together with the corresponding validation performance.

Whole Slide Image Maps

Tile-level segmentation masks were aggregated to produce an axon density map (panel B), where pixel intensity reflects the fraction of foreground axon pixels per spatial bin. We chose a spatial bin of 70 × 70 μm to maintain anatomical detail while keeping the image size manageable. The spatial bin size can be easily adjusted in the software application. Axon orientation was quantified using a grid-based structural-tensor approach applied to the binary whole-slide axon segmentation mask. The mask was subdivided into a regular grid of non-overlapping square cells (again with a spatial bin of 70 × 70 μm), and for each cell the following procedure was performed. Image gradients were computed using the Scharr operator (Scharr, 2007 ). The outer-product structure tensor was then assembled per pixel as:

and each component was smoothed with a Gaussian kernel of standard deviation σ = 2 pixels to integrate local neighborhood information. The cell-level tensor was obtained by averaging over all axon-positive pixels within the cell. Eigenanalysis of yielded the dominant and secondary eigenvalues, with corresponding eigenvectors defining the primary orientation and perpendicular secondary orientation , both expressed in the range [0°, 180) via

Alignment strength per cell was quantified as the 2D fractional anisotropy (FA)-analog coherency:

ranging from 0 (isotropic, no preferred direction) to 1 (perfectly aligned). Final output maps consisted of an axon density map, a coherence map and an HSV-encoded summary image in which hue encodes , saturation encodes coherency, and value indicates axon presence (set to 1 if axons are present within the grid cell and to 0 in empty cells).

Implementation

The segmentation pipeline was implemented using PyTorch 2.x (Paszke et al., 2019 ) with PyTorch Lightning for training orchestration. Data augmentation utilized the Albumentations library (Buslaev et al., 2020 ), and experiment tracking was performed using Weights & Biases ( n.d. ). Training was conducted on an NVIDIA RTX 2000 Ada GPU (8 GB VRAM) with batch size 16. Reproducibility was ensured by setting a global random seed (42) and enabling deterministic algorithms. Model selection was based on minimum validation loss, with the top three models and the most recent checkpoint preserved.

Baseline Comparison (nnU-Net)

As an external baseline we trained nnU-Net v2 (Isensee et al., 2021 ), a widely used self-configuring segmentation framework, on the identical data. To keep the comparison controlled, nnU-Net received the same 256 × 256 patches and the same data partition as our model: the 17,328 un-augmented training patches and 2,032 validation patches were supplied with our exact train/validation split imposed as a custom single fold, and the 2,032 test patches were held out. The 2D configuration was trained with nnU-Net’s default self-configuring pipeline (automated planner, 256 × 256 patch size matching our field of view, 1000 epochs, the framework’s own on-the-fly data augmentation), and its standard best-validation checkpoint (selected on validation pseudo-Dice) was used for inference. No test-time post-processing was applied, matching our model. nnU-Net predictions were binarized at 0.5 and evaluated with the identical stratified per-tile protocol (Sect. 2.7 ), so the resulting metrics are directly comparable to those of our Attention U-Net. For the qualitative whole-tile comparison, each model was additionally applied to the full 512 × 512 test tiles (nnU-Net through its standard sliding-window inference).

Results

Training Performance

The full model (7.9 million trainable parameters) was trained on a total of 86,330 augmented patches (5,415 gradient steps per epoch, batch size 16) for a maximum of 100 epochs with early stopping monitored on validation loss (patience: 5 epochs, minimum improvement threshold: ).

Training converged rapidly. The Focal Tversky Loss decreased steeply during the first two epochs and reached its minimum on the validation set at epoch 2 (val/loss 0.0528 ; Fig. 1 ). Epochs 3 through 7 showed no further improvement in validation loss, triggering early stopping after a total of 8 epochs (ca. 5.5 h on an NVIDIA RTX 2000 Ada). The best checkpoint (epoch 2) was restored for all downstream evaluations.

All validation metrics (Precision, Recall, F1, IoU) tracked the validation loss improvement closely, with their best values reached at or near the same checkpoint (Fig. 2 ). The recall ( at best checkpoint) consistently exceeded precision ( ), reflecting the asymmetric Focal Tversky loss design which penalizes false negatives more heavily than false positives ( ).

Ablation Study: Loss Weighting and Attention Gates

To empirically isolate the two important design choices — the Focal Tversky weighting and the attention-gated skip connections — we performed two ablations under the same held-out test set and stratified per-tile evaluation used above (8 tiles, threshold 0.5; Supplementary Results 1 ). All other settings (architecture, optimizer, augmentation, and the focusing parameter γ = 2.0) were held fixed.

Loss Weighting

We swept the Tversky weightings along α + β = 1 from β = 0.1 to β = 0.9 (Supplementary Tables 1 and Supplementary Fig. 3 ). Test Dice peaked at β = 0.6 (Dice = 0.715), the value used in the main model, but the top of the curve was a broad plateau: settings from β = 0.4 to β = 0.8 all scored within approximately 0.01 Dice of the maximum, a difference smaller than the tile-to-tile standard deviation. As expected, precision decreased and recall increased monotonically with β. We therefore regard β = 0.6 as the highest-scoring setting; our selection was additionally motivated by a deliberate, task-motivated preference for recall over precision, since for axon-density quantification a missed axon was considered more detrimental than a slight over-segmentation.

Attention Gates

At the selected loss configuration (α = 0.4, β = 0.6) we compared the network with and without the attention-gated skip connections, training three independent seeds of each (Supplementary Fig. 4 ). Enabling the attention gates improved mean test Dice from 0.699 to 0.715 (mean ΔDice = + 0.017), and the improvement was consistent across all three seeds and all five metrics (Dice, IoU, precision, recall, F1), at a cost of only about 0.09 M additional parameters (approximately 1% of the model). The effect is thus small but consistent and essentially free in model size, and we observed no downside to including the gates.

Test Set Performance

On the held-out test set the model achieved a mean pixel-wise F1/Dice of and an IoU of (mean SD across the 8 test tiles; Table 2 ). Mean recall ( ) exceeded mean precision ( , consistent with the design intent of the Focal Tversky Loss to prioritize recall and minimize missed axon pixels (Fig. 3 ).

Performance varied considerably across individual test tiles. The highest Dice score was achieved on tile 3_4 (Dice 1.00 , IoU 1.00 ), while the lowest was observed on tile 3_14 (Dice 0.392 , IoU 0.244 ). Note that tile 3_4 was a background tile, explaining the high evaluation performance. The low performance on tile 3_14 can be explained by the very low staining quality, making axon detection difficult also for human raters.

Notably, the overall performance of 0.717 Dice is encouraging in the context of the substantial inter-rater variability inherent in manual axon annotation (Table 3 ; Fig. 4 ), suggesting that the model approaches the practical upper bound of reproducible annotation quality.

Inter-Rater Agreement

F1-scores between the two raters were 0.637 for tile 1_3 and 0.671 for tile 4_3. The visual differences between the two rater masks (Fig. 4 ) further illustrate the reasons: both raters segment similarly the major axons, while thin branches and low-contrast fibers are demonstrating differences. Another valid point of difference is each rater’s perception of axonal width.

Critically, the model’s agreement with each individual rater (F1: 0.736 for rater 1 and 0.703 for rater 2 on tile 1_3; 0.723 for rater 1 and 0.681 for rater 2 on for tile 4_3) met or exceeded the rater-to-rater agreement on both tiles. These results demonstrate that the model’s predictions are as consistent with each human annotation as the two human annotators are with each other, confirming that the model operates within the intrinsic variability of the annotation task.

Leave-One-WSI-Out Cross-Validation

LOWO models were evaluated on checkpoints at validation loss minima, which occurred between epoch 1 and epoch 10 across folds (median: epoch 4–5), again demonstrating rapid convergence. Across all 12 folds, the LOWO analysis yielded a mean Dice of (Table 4 ; Fig. 5 ), with limited inter-fold variability (range: 0.687–0.725). The low standard deviation suggests that no single training WSI is disproportionately critical to performance. The highest fold performance was achieved when WSI 2_3 was excluded (Dice 0.725 ), while the lowest was observed when WSI 3_15 was left out (Dice 0.687 ).

Qualitative Results

Figure 6 shows 7 test set tiles (the background tile 3_4 with Dice 1.000 is not shown) alongside their corresponding model-predicted axon segmentation masks. The model correctly recovered the densely packed axon network in white matter regions and successfully suppressed non-axonal structures such as cell bodies and background artefact. Prediction accuracy was reduced in regions with increased background staining and reduced axon density.

Figure 7 shows a color-coded prediction overlay for tile 4_1, which achieved Dice 0.773 on the test set. The overlay depicts true-positive (green), false-negatives (red), and false-positives (orange) detections, enabling direct visualization of prediction errors within the tissue context. While axons are delineated with high accuracy in general, fine dendrite structures are often false-negative. False-positives arise mainly in regions with ambiguous stain or due to background noise.

To demonstrate applicability beyond isolated tiles, we applied the trained model to a complete Bielschowsky-stained whole-slide image (WSI 1_3) and derived three spatially resolved maps, downsampled to 70 μm resolution (Fig. 8 ). The resulting map captures the expected regional anatomy: white matter regions exhibit uniformly high axon density, whereas the cortical gray matter and the MS lesion area show reduced density, consistent with the known cytoarchitecture of the tissue.

The orientation map (panel D) reveals varying fiber direction domains that correspond to white matter tracts. The coherency map (panel C) quantifies the degree of directional alignment of local axonal architecture. Highly organized fiber bundles correspond to brighter regions, whereas crossing-fiber regions, cortex, and lesion areas exhibit lower coherency (darker regions).

Baseline Comparison (nnU-Net)

Across all eight test tiles our Attention U-Net attained a higher mean Dice than nnU-Net (0.717 ± 0.162 vs. 0.555 ± 0.245; IoU 0.584 ± 0.201 vs. 0.417 ± 0.200; paired Wilcoxon signed-rank on per-tile Dice, p = 0.023, n = 8; Supplementary Results 2 , Supplementary Table 2 ). This full-set gap was strongly influenced by the background tile 3_4, whose empty ground truth makes the pixel Dice degenerate: our model predicted no axons there (Dice 1.000), whereas a few false-positive pixels drove nnU-Net’s patch-pooled score to 0.000, even though its whole-tile prediction on this tile is likewise essentially empty (Supplementary Fig. 5 ). To report the comparison fairly, we therefore also evaluated the seven axon-bearing tiles, excluding this degenerate background tile. On these tiles the two models were closer, but our model remained ahead on Dice (0.677 ± 0.130 vs. 0.634 ± 0.135; p = 0.047, n = 7), with a complementary error profile: the Attention U-Net achieved higher recall (0.719 vs. 0.606) whereas nnU-Net was more precise (0.749 vs. 0.676). Per-tile metrics for both models are given in Supplementary Tables 3 , and a whole-tile qualitative comparison in Supplementary Fig. 5 . Notably, our model reached this performance with far fewer parameters (7.9 M vs. approximately 48 M for nnU-Net’s self-configured network).

Discussion

We present an attention U-Net trained to segment axons in Bielschowsky silver-impregnated histological sections. The achieved mean pixel-wise F1/Dice score and Leave-one-WSI-out cross-validation on the held-out test set, spanning four brains, multiple tissue classes, and both normal and demyelinated tissue, demonstrated generalizable and robust model performance. Furthermore, we demonstrate that the binary segmentation output can be aggregated to full WSI-level maps of axon density, dominant fiber orientation, and local fiber coherency (Fig. 8 ), establishing a practical pipeline from raw silver-stained scan to spatially resolved quantitative readouts.

For thin, densely packed structures such as axons in silver-stained tissue, human-level agreement is itself limited. In our own annotation experiment on two representative tiles, rater-to-rater F1 was 0.637 and 0.671, whereas the model-to-rater agreements were 0.681–0.736 (Table 3 ), consistently at or above the human ceiling. Comparisons for similar tasks are rare in the literature. Golabchi et al. ( 2010 ) reported stereological inter-rater variability for axon counting in gray matter. Deng et al. ( 2021 ) reported a Dice coefficient of 0.78 between two raters for pixel-wise axon segmentation in paraphenylenediamine-stained optic nerve sections, and Kumar et al. ( 2020 ) reported an average aggregated Jaccard index of 0.653 in pixel-wise nuclei segmentation in hematoxylin-eosin sections. Our model showed a Dice score of 0.72, indicating that it operates within the range of inter-rater disagreement and suggesting that the principal remaining source of error is annotation ambiguity rather than model-induced limitations.

Two design choices were central to the model’s performance: attention-gated skip connections and the Focal Tversky Loss. The attention mechanism allows the decoder to selectively weight encoder features according to relevance at each spatial location, suppressing responses in background regions while preserving detail in axon-dense areas. This is particularly valuable in Bielschowsky-stained tissue, where heterogeneous background silver deposition can produce false-positive activation in a standard U-Net. The Focal Tversky Loss addresses the severe class imbalance inherent in axon segmentation (approximately 30% foreground pixels) by penalizing false negatives more heavily than false positives ( ) and applying a quadratic focusing factor ( ) that concentrates learning on difficult examples. The consequence is visible in the test set metrics: mean recall ( 0.754 ) consistently exceeded mean precision ( 0.717 ), indicating that the model errs toward detecting axons rather than missing them, which we consider a less detrimental error pattern for density quantification, where missed axons directly reduce density estimates. To empirically substantiate these two choices, we performed a Tversky-weighting sweep and an attention on/off comparison (Sect. 3.2 ). The sweep identified β = 0.6 as the highest-scoring weighting, situated at the top of a broad plateau, and our choice additionally reflects a task-motivated preference for recall over precision. The attention gates provided a small but consistent improvement across seeds at negligible parameter cost. We therefore regard both as sensible, empirically supported design choices.

Beyond tile-level evaluation, the WSI-level density, orientation, and coherency maps derived from the model output (Fig. 8 ) illustrate the practical utility of the pipeline for research and lab-routine applications. The density map recapitulates the expected gross anatomy of the tissue block: compact white matter with uniformly high axon density, a cortical gray matter ribbon with lower and more variable density, and visibly reduced density in MS-lesion areas. The structural tensor–derived orientation and coherency maps provide an orthogonal readout, where fiber coherency is increased in organized, highly-directional WM tracts and decreased at GM boundaries and lesion borders, consistent with the known disruption of myeloarchitecture in those regions.

Our work comes with some limitations. Despite training and test tiles were resected from different brains, different anatomical locations and were prepared at different time points by more than one technician, they originate from a single disease context (multiple sclerosis), and a single staining protocol was followed. This might constitute a limiting factor on generalization over different silver impregnation protocols. Expanding the training dataset to include additional brains, other disease contexts (e.g., Alzheimer’s disease, traumatic brain injury), and sections stained in different laboratories would improve generalization and allow a more rigorous evaluation of domain robustness.

Regarding the practical scope of the model, several boundaries should be kept in mind. All training data were digitized at a single resolution (0.273 μm per pixel); when we applied the trained model, without retraining, to Bielschowsky sections digitized at a coarser resolutions (0.35 μm per pixel and 0.70 μm per pixel), the segmentations remained visually convincing (Supplementary Results 3 , Supplementary Fig. 6 ). We did not, however, annotate data at other resolutions and therefore cannot quantify performance there, and behavior under more pronounced changes in resolution is unknown. Because our training tiles deliberately include sections prepared by different laboratory technicians at different time points, and because we apply strong data augmentation, we expect the model to transfer reasonably to Bielschowsky material from other laboratories acquired at a comparable resolution; we have, however, no data from other laboratories to confirm this, and generalization to other silver-impregnation stains was not evaluated. Where such domain shifts prove too large, a foundation-model-based approach (e.g., fine-tuning a large pretrained segmentation model) is a promising future direction for broadening applicability.

The annotation set is relatively small ( n=33 tiles), which might constrain both the diversity of the training set and the statistical power of the evaluation. Additionally, the complexity of the annotation task (particularly for thin dendrite branches and regions of non-specific background staining) propagates directly into the reported performance metrics. Involvement of more primary annotators for the segmentations, would potentially provide a better-grounded ceiling for model performance and help disentangle model error from annotation noise. Because dense axon annotation is so labor-intensive, each tile was segmented only once, by a single rater (except the two tiles of the dedicated inter-rater experiment); we therefore lack repeated independent annotations of the modelling tiles, and the inter-rater ceiling was estimated from two tiles rather than the full set.

Despite strong mean performance, it is important to note that test Dice scores ranged widely from ~ 0.4 to 1 across the eight test tiles, reflecting the influence of staining quality and regional complexity on segmentation performance. Uniformly stained tissue regions with clearly distinguishable axonal structures achieved the highest performance metrics, whereas tiles with severe background stain, complex structural details and lacking contrast exhibited substantially lower segmentation accuracy. However, similar challenges also affect manual segmentation, where poor staining quality (e.g., non-specific silver precipitates, stain inhomogeneities, and intricate tissue morphology), likewise complicate accurate annotation. To our knowledge, no axon segmentation tool tailored to Bielschowsky-stained tissue sections exists, so a task-specific benchmark was not available; we therefore compared our model against nnU-Net v2, a strong self-configuring segmentation framework, trained on the identical data (Sect. 3.7 ).

On this benchmark our Attention U-Net was on Dice modestly above the nnU-Net across both the full eight-tile test set and the seven axon-bearing test tiles, while using far fewer parameters (7.9 M versus approximately 48 M). The clearest advantage of our model was on the white-matter tile 2_3 (Dice 0.805 versus 0.630), where axons run largely perpendicular to the sectioning plane and therefore appear as small cross-sectional profiles rather than elongated fibers. Given the small number of test tiles ( n = 8), this comparison should be read as an honest benchmark rather than a definitive ranking, and it indicates that our compact, task-specific model performs on par with, or slightly better than, a fully self-configuring pipeline at a fraction of the model size.

Methodologically, we acknowledge that we trained the model in a small dataset, which bears the intrinsic risk that its performance depends on individual training examples. However, the LOWO analysis indicates that model performance is not strongly dependent on any single training WSI, as shown by the narrow Dice range across folds (0.687–0.725, SD 0.010). This suggests that the small training dataset provides complementary coverage of relevant tissue variation, and that robust segmentation can be achieved from a limited number of annotated tiles when extensive patch extraction and data augmentation are used.

Summary

In this study, we present a deep learning-based tool for automatic axon segmentation in Bielschowsky-stained tissue sections, achieving a pixel-wise F1/Dice of 0.717 ± 0.162 on a held-out test set spanning four brains, multiple tissue classes, and both normal and demyelinated tissue — performance that meets or exceeds human inter-rater agreement. The attention U-Net architecture, trained with Focal Tversky Loss to address class imbalance, produces tile-level and whole-slide segmentation maps from which spatially resolved axon density, dominant fiber orientation, and local fiber coherency metrics can be derived. Leave-one-WSI-out analysis confirmed robustness to individual training examples (Dice range: 0.687–0.725), supporting practical utility even with limited annotated data. The pipeline is released as an open-source software package for straightforward deployment on tiled images and whole-slide scans. Future integration with multimodal registration frameworks would enable direct spatial alignment between the derived WSI orientation maps and ex-vivo or in-vivo MRI datasets, establishing a quantitative histological ground truth for microstructure imaging models and opening new avenues for bridging mesoscopic and macroscopic analyses of axonal pathology.

Supplementary Information

Below is the link to the electronic supplementary material.

Supplementary Material 1 (DOCX 97.0 MB)

Acknowledgements

The authors thank Federico Spagnolo for valuable insights and helpful input throughout the development of this work. His feedback and constructive suggestions improved the quality and clarity of the manuscript.

Author Contributions

Lukas Schönenberger: Conceptualization, Methodology, Software, Formal analysis, Writing – original draft. Laurin Egli: Data curation, Writing – review & editing. Dimitrios Gkotsoulias: Methodology, Writing – review & editing. Christine Stadelmann: Resources, Writing – review & editing. Cristina Granziera: Resources, Writing – review & editing, Supervision, Funding acquisition.

Funding

Open access funding provided by University of Basel

Data Availability

Source code (version 1.0.1) is archived on Zenodo and available at https://doi.org/10.5281/zenodo.19455087. The trained model is available on Hugging Face at https://huggingface.co/lukas-schoenenberger/BasNet/tree/main. Annotated training data are deposited on a private Hugging Face repository at https://huggingface.co/datasets/lukas-schoenenberger/Basnet-dataset and are available from the corresponding author upon reasonable request with the permission of the German MS Brain Bank of the Competence Network Multiple Sclerosis.

Informed Consent

Not applicable. The study used human post-mortem brain tissue obtained through the German MS Brain Bank.

Research Involving Human Participants and/or Animals

This study used human post-mortem brain tissue provided by the German MS Brain Bank of the Competence Network Multiple Sclerosis. The study was approved by the ethical review committee of the University Medical Center Göttingen.

Competing interests

The authors declare no competing interests.

Footnotes

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Christine Stadelmann and Cristina Granziera contributed equally.

Contributor Information

Lukas Schönenberger, Email: lu.schoenenberger@unibas.ch.

Cristina Granziera, Email: cristina.granziera@unibas.ch.

References

Abdollahzadeh, A., Belevich, I., Jokitalo, E., Sierra, A., & Tohka, J. (2021). DeepACSON automated segmentation of white matter in 3D electron microscopy. Communications Biology, 4 (1), Article 179. 10.1038/s42003-021-01699-w

Abraham, N., & Khan, N. M. (2018). A novel focal Tversky loss function with improved Attention U-Net for lesion segmentation. arXiv . 10.48550/arXiv.1810.07842

Bielschowsky, M. (1902). Die Silberimpragnation der Achsenzylinder. Neurologisches Zentralblatt , 21 , 579–584.

Buslaev, A., Parinov, A., Khvedchenya, E., Iglovikov, V. I., Druzhinin, M., & Kalinin, A. A. (2020). Albumentations: Fast and flexible image augmentations. Information , 11 (2), 125. 10.3390/info11020125

Deng, W., Hedberg-Buenz, A., Soukup, D. A., et al. (2021). AxonDeep: Automated optic nerve axon segmentation in mice with deep learning. Translational Vision Science & Technology, 10 (14), Article 22. 10.1167/tvst.10.14.22

Fathy, Y. Y., Jonkman, L. E., Bol, J. J., et al. (2022). Axonal degeneration in the anterior insular cortex is associated with Alzheimer’s co-pathology in Parkinson’s disease and dementia with Lewy bodies. Translational Neurodegeneration, 11 (1), Article 52. 10.1186/s40035-022-00325-x

Gkotsoulias, D. G., Schonenberger, L., Leupold, J., et al. (2025). Bridging mesoscopic and microscopic scales in multiple sclerosis: Post mortem brain block multicontrast 9.4T MRI and histology quantification. STAR Protocols, 6 (4), Article 104203. 10.1016/j.xpro.2025.104203

Golabchi, F. N., Brooks, D. H., Hoge, W. S., De Girolami, U., & Maier, S. E. (2010). Pixel-based comparison of spinal cord MR diffusion anisotropy with axon packing parameters. Magnetic Resonance in Medicine, 63 (6), 1510–1519. 10.1002/mrm.22337

Isensee, F., Jaeger, P. F., Kohl, S. A. A., Petersen, J., & Maier-Hein, K. H. (2021). nnU-Net: A self-configuring method for deep learning-based biomedical image segmentation. Nature Methods , 18 (2), 203–211. 10.1038/s41592-020-01008-z

Kevenaar, J. T., & Hoogenraad, C. C. (2015). The axonal cytoskeleton: From organization to function. Frontiers in Molecular Neuroscience, 8 , Article 44. 10.3389/fnmol.2015.00044

Kumar, N., Verma, R., Anand, D., et al. (2020). A multi-organ nucleus segmentation challenge. IEEE Transactions on Medical Imaging, 39 (5), 1380–1391. 10.1109/TMI.2019.2947628

Loshchilov, I., & Hutter, F. (2017). Decoupled weight decay regularization. arXiv . 10.48550/arXiv.1711.05101

Nestor, S. L. (2008). Techniques in neuropathology. In J. D. Bancroft & M. Gamble (Eds.), Theory and practice of histological techniques (6th ed., pp. 365–403). Churchill Livingstone Elsevier. 10.1016/B978-0-443-10279-0.50026-9

Oktay, O., Schlemper, J., Le Folgoc, L. (2018). Attention U-Net: Learning where to look for the pancreas. arXiv . http://arxiv.org/abs/1804.03999

Paszke, A., Gross, S., Massa, F. (2019). PyTorch: An imperative style, high-performance deep learning library. arXiv . http://arxiv.org/abs/1912.01703

Ronneberger, O., Fischer, P., & Brox, T. (2015). U-Net: Convolutional networks for biomedical image segmentation. arXiv . 10.48550/arXiv.1505.04597

Scharr, H. (2007). Optimal filters for extended optical flow. In B. Jahne, R. Mester, E. Barth, & H. Scharr (Eds.), Complex motion (Vol. Vol. 3417, pp. 14–29). Springer. 10.1007/978-3-540-69866-1_2

Stadelmann, C., Timmler, S., Barrantes-Freer, A., & Simons, M. (2019). Myelin in the central nervous system: Structure, function, and pathology. Physiological Reviews , 99 (3), 1381–1431. 10.1152/physrev.00031.2018

Tversky, A. (1977). Features of similarity. Psychological Review , 84 (4), 327–352. 10.1037/0033-295X.84.4.327

Weights & Biases. (n.d.). Weights & Biases: The AI developer platform. Retrieved March 16, 2026 from https://wandb.ai/site/

Zaimi, A., Wabartha, M., Herman, V., Antonsanti, P. L., Perone, C. S., & Cohen-Adad, J. (2018). AxonDeepSeg: Automatic axon and myelin segmentation from microscopy data using convolutional neural networks. Scientific Reports, 8 (1), Article 3816. 10.1038/s41598-018-22181-4

Zuiderveld, K. J. (1994). Contrast limited adaptive histogram equalization. In P. S. Heckbert (Ed.), Graphics gems (pp. 474–485). Academic Press. 10.1016/B978-0-12-336156-1.50061-6

Associated Data

Supplementary Materials

Supplementary Material 1 (DOCX 97.0 MB)

Data Availability Statement

Source code (version 1.0.1) is archived on Zenodo and available at https://doi.org/10.5281/zenodo.19455087. The trained model is available on Hugging Face at https://huggingface.co/lukas-schoenenberger/BasNet/tree/main. Annotated training data are deposited on a private Hugging Face repository at https://huggingface.co/datasets/lukas-schoenenberger/Basnet-dataset and are available from the corresponding author upon reasonable request with the permission of the German MS Brain Bank of the Competence Network Multiple Sclerosis.