Deep Learning-Based Quantification of Pulmonary Hemosiderophages in Cytology Slides

Exercise-induced pulmonary hemorrhage (EIPH) is a common condition in sport horses with negative impact on performance. Cytology of bronchoalveolar lavage fluid by use of a scoring system is considered the most sensitive diagnostic method. Macrophages are classified depending on the degree of cytoplasmic hemosiderin content. The current gold standard is manual grading, which is however monotonous and time-consuming. We evaluated state-of-the-art deep learning-based methods for single cell macrophage classification and compared them against the performance of nine cytology experts and evaluated inter- and intra-observer variability. Additionally, we evaluated object detection methods on a novel data set of 17 completely annotated cytology whole slide images (WSI) containing 78,047 hemosiderophages. Our deep learning-based approach reached a concordance of 0.85, partially exceeding human expert concordance (0.68 to 0.86, mean of 0.73, SD of 0.04). Intra-observer variability was high (0.68 to 0.88) and inter-observer concordance was moderate (Fleiss’ kappa = 0.67). Our object detection approach has a mean average precision of 0.66 over the five classes from the whole slide gigapixel image and a computation time of below two minutes. To mitigate the high inter- and intra-rater variability, we propose our automated object detection pipeline, enabling accurate, reproducible and quick EIPH scoring in WSI.


Related Work
To date, the topic of hemosiderophage classification and quantification has not been approached using computer vision methods. However, there have been numerous studies in the past decades with the goal of detecting cells, nuclei and mitotic figures for multiple modalities like digital fluorescence microscopy and histopathology [17][18][19] . Historically, this started as hand-crafted low-level feature extraction [20][21][22] . With the recent advent of deep learning-based techniques 23 these methods transitioned into modern end-to-end optimised object detection algorithms like Faster-RCNN 24 , SSD 25 or RetinaNet 26 . Their underlying end-to-end optimisation approach is the foundation of their success in object detection challenges for natural images like PASCAL VOC 27 and MS COCO 28 where no classical approach could outperform a modern deep learning-based object detection method 29 since 2014. The aim of object detection algorithms is to predict the bounding box as well as a class for multiple objects irrespective of the scale or a partial occlusion of the objects. These methods have generated state of the art results in the fields of pedestrian-, face-and car-detection and are used in state of the art autonomous vehicles as well as the interpretation of satellite images 27,28,30 . Regarding the field of digital pathology object detection, the review by Litjens et al. 31 reveals that no one had implemented deep learning-based object detection methods for the evaluation of medical images as of 2017. In contrast, they mention that sliding window approaches in combination with a deep learning-based classification network or U-Net-like segmentation architectures 32 are being commonly used. The frequent use of U-Net in particular is quite remarkable since segmentation provides no means of separating touching or overlapping objects and these methods highly rely on post-processing steps for the task of separation. Additionally, in the case of U-Net, the architectures are computationally more complex due to their encoder-decoder architecture. Moreover, these networks require a pixel-wise annotation mask for obtaining better results, which is time-consuming compared to the relatively simple and fast creation of bounding box annotations needed for object detection methods. Ferlaino et al. 33 used deep learning-based object detection on fully annotated multiclass WSI. For this, they employed RetinaNet 26 for nuclei detection and a separate, not end-to-end trainable, network for nuclei classification.
Modern object detection approaches can be categorised into the two major categories of single stage and two stage algorithms. In single stage setups, the task of detection and classification is solved in one single network, examples are YOLO 34 , SSD 25 or RetinaNet 26 . In two stage algorithms, the task of detection is solved by the use of a region proposal network (RPN) 24 in the first stage and then classified using an additional network in a subsequent stage. While two stage detection is more accurate in general, the single stage methods yield the better ratio of accuracy and inference speed 35 . This trade-off between speed and accuracy is crucial when analysing WSI with billions of pixels. In this work, we used RetinaNet as a starting point for analysing EIPH on WSI because its architecture is straightforward, easy to modify and adapt for WSI analysis.

Material
Our research group built a data set of 17 cytological slides of equine bronchoalveolar lavage fluid. The slides were prepared by cytocentrifugation and stained for iron content with Prussian Blue (n = 10) or Turnbull's Blue (n = 7) which results in an identical colour pattern. Digitalisation of the glass slide was performed using a linear scanner (Aperio ScanScope CS2, Leica Biosystems, Germany) at a magnification of 400× (resolution: . ). Finally, all macrophages on each slide were annotated by a veterinary pathologist. All bronchoalveolar lavage fluids were obtained from horses with clinical signs of lower respiratory tract disease during routine diagnostic service for therapeutic reasons. Written informed consent was obtained from the owners. Therefore, no animal was harmed for the construction of this data set. Individual case histories were not considered in the present study and all data we received was anonymised by the routine diagnostic service in advance. Using the open source software solution SlideRunner 36 , we were able to build a database that includes the annotations for each hemosiderophage on the slides with their corresponding grade. This was done by first annotating all pulmonary macrophages and afterwards classifying them into their corresponding grade. The scoring system for hemosiderophages was introduced by Golde et al. 4 and consists of five classes: It ranges from zero (no intracytoplasmic blue coloured pigment) to four (cell filled with hemosiderin; dark blue throughout cytoplasm). The final score was calculated by the www.nature.com/scientificreports www.nature.com/scientificreports/ method of M. Y. Doucet and L. Viel 16 which is an adaptation of Golde et al. 4 to be used for horses. In this scoring system, three hundred alveolar macrophages were first graded from zero to four, then the total number per grade was divided by three and multiplied with the corresponding grade. The resulting total hemosiderin score (THS) thus ranges from zero to four hundred. If the score is higher than 75 then the diagnosis pulmonary hemorrhage is considered to be confirmed. The completely annotated data set consists of 17 slides and covers an area of 1,266 mm 2 (mean = 74 mm 2 , SD = 9 mm 2 ) containing 78,047 labelled cells (mean = 4,591, SD = 3,389) (see Table 1) making it the largest published data set of hemosiderophages and one of the largest of WSI. This novel data set allows us to perform object detection on whole hemosiderophages slides for the first time.

Methods
The research was carried out in accordance with the Code of Ethics of the World Medical Association (Declaration of Helsinki) and the guidelines of the institutions conducting the experiments.
The aim of this work was to develop and compare algorithmic approaches for predicting the hemosiderophages score of WSI. In order to assess how challenging the classification of single hemosiderophages (CoSH) is, we investigated two methods considering the single cell labels as a classification and as a regression task. We then compared the results with human performance. Additionally, we present methods for multi-class WSI analysis (MCWSA). Here, we adopted state of the art deep learning-based object detection and regression approaches. We used a support vector regression to draw a baseline. To compensate for the sparse cell distribution, we introduce a novel quadtree-based sampling approach to train the object detection networks.
Human performance evaluation. In order to compare our algorithmic approaches with human recognition performance, we investigated the accuracy and reproducibility of nine cytology experts in labelling single cell pulmonary hemosiderophages. We divided them into three groups according to their qualification and experience with BAL cytology. Each group contained three participants: • (E)xpert: Veterinary pathologists or clinician with high degree of experience in BAL cytology.
• (P)rofessional: Professional clinician or pathologist with basic experience in BAL cytology.
• (B)eginner: General skills in cytology, but no experience with BAL cytology in particular.
To evaluate the human inter-and intra-observer variability for single cell classification, we extracted two test sets containing 1,000 cells each. For test set 1, the images were randomly selected among the labelled cells resulting in a representative distribution. Test set 2 contained 1,000 cells with a balanced distribution of 200 cells per grade.
Each of the nine cytology experts was asked to classify two thousand cells from the single cell test set 1 and 2. We did not set a time limit to perform this task. In order to measure the intra-observer variability, they were asked to classify all cells again two weeks after the initial assessment. The participants were instructed to perform classification according to the methods published by Doucet et al. 16 .
Sampling strategy. Taking into account that not all slides contain hemosiderophages of grade three and four, we used the same fourteen slides to train and validate. However, we used the upper half of each image for training and the lower half for validation in order to prevent over-fitting. Three separate slides were selected as hold out test set slides.  Table 1, cell grade 3 and 4 rarely occur on some of the WSI. For example, slide 14 includes only one grade 4 and eight grade 3 hemosiderophages. This means that with an image size of 35,999 × 34,118 pixels and random sampling with a patch size of 1024 × 1024 pixels, the chance to sample the grade 4 cell is only 0.08% percent.
Two-stage cluster sampling strategies. For this sampling strategy, we clustered all cells from one WSI on the basis of their grade. For training, we randomly selected one of those clusters and chose one of the cells within that cluster by chance. Then a patch is randomly shifted in the direct proximity of that cell and the area is sampled for training.
Generic quadtree sampling strategies. We developed a novel sampling strategy for microscopy images based on a quadtree in order to consider the probability of occurrence of cells as well as their neighbouring cells (see Fig. 1 center). At each level of the quadtree (depth of the tree can be customised), we saved the cells, their corresponding sampling probability and their grade. As seen in Fig. 1 (center), at each level of the quadtree, we have up to four nodes. One constraint for the tree while it is being created was that there must be at least three hundred cells in each node. One other option would be that the size of the final node must be identical to the training patch size (e.g. 1024 × 1024 pixels). In contrast to the sampling strategy described in the previous section, we can sample at nodes without any cells by defining a minimum probability. To train our networks, we created a quadtree with a depth of three. To create a training sample, we randomly traversed the quadtree according to the sampling probability of the cells. Figure 1 visualises this novel sampling approach. At the first level, the image is divided into four nodes with the sampling probabilities of 35.3%, 32.4%, 13.9% and 18.4% (clockwise). In this example, the top right node was selected by chance and was traversed further. This process was repeated until the final node at level three was reached and one patch was extracted for training.

Single cell classification (CoCH).
The hemosiderophages score is based on a subjective, semi-quantitative method in which each cell in a selected region of the WSI is assigned one out of five grades (ranging from zero to four). However, this quantised grading system does not reflect the biological nature since there is a continuous gain of iron in the hemophages as opposed to a stepwise rise. To take this continuous increase into account, we propose a regression-based cell score estimation. We then compare the result to the classification approach mimicking the human scoring system.
Classification. For the cell-based classification task, we used a compact ResNet-18 Architecture 37 pre-trained on ImageNet 38 with a fully connected two layer classification head and a final softmax activation. The cells used for training and validation were extracted according to the proposed quadtree-based sampling strategy. The Network was trained in two stages with the Adam optimiser and a maximal learning rate schedule of 0.01. Categorical www.nature.com/scientificreports www.nature.com/scientificreports/ cross entropy was used as the loss function. First, we trained only the classification head for three epochs, afterwards we fine-tuned the complete network for an additional twenty epochs until convergence was reached.
Regression. As stated, the hemosiderin absorption is a continuous process which is mapped to a discrete grading system. To take his continuity into account, we developed a network with a regression head and a final scaled sigmoid activation which predicts continuous values in a range of −0.5 to 4.5. This compensates for the implementation instability for sigmoid activations close to zero and one. The main focus of the experiment was to estimate the intra-grade confusion and increase the human interpretability of the results. This modification enables the network to predict decimal values between any two grades given that the cell has features supporting two grades, which is not possible with a classification approach (see Fig. 2). The network and training schedule were applied as described in the single cell classification paragraph. The mean squared error was used as the loss function.
Object detection-based WSI score estimation (MCWSA). Besides investigating pure classification performance on single cells where the coordinates are previously known, the actual task in diagnostics is the estimation of scores on complete WSIs or subparts thereof. Object detection networks mimic human expert behaviour by both detecting and classifying the cells and calculating the score afterwards. One object detection approach with a good accuracy-speed trade-off is RetinaNet 26 which is a single, unified network composed of a backbone network for feature extraction (see Fig. 3a). A feature pyramid network (FPN) 39 is built on top of the feature extractor to generate rich, multi-scale features by combining low-resolution with semantically strong features and high-resolution with semantically weak features (see Fig. 3c). On each layer of the FPN, a classification subnet and a regression subnet are called to make predictions (see Fig. 3d,e). The classification head predicts the probability of the target object's presence at each spatial position for each anchor. Anchors are defined by the scale and aspect ratio to match the targeted objects on each spatial position. To compensate for the class imbalance,   Figure 2. Cell-based regression results on the test data set visualised as a density histogram for the predicted scores. As an example, both cells in the middle are labelled with grade two and the regression model assigned very different scores to both, which is also clearly comprehensible from the visual appearance of the cell. www.nature.com/scientificreports www.nature.com/scientificreports/ focal loss 26 was employed during training. The bounding box regression subnet (see Fig. 3f) is generally built in a similar fashion as the classification head but was trained with smooth L1 loss and predicted four coordinates (x-offset, y-offset, width, height) for each box if a corresponding anchor box existed.
We have modified the RetinaNet architecture in three significant ways to further optimise it for hemosiderophage WSI analysis. Firstly, we added an additional regression head which predicts the hemosiderophages score for each hemosiderophage (see Fig. 3f). This had the intent to increase the human interpretability of the results. As the loss function for the cell-based regression head, mean squared error was used. Secondly, to utilise the features extracted from the RetinaNet backbone, we fitted an additional regression head on top of the ResNet-18 feature extractor for patch-wise hemosiderophages score prediction. This process is further described in the later section deep learning-based regression and visualised in Fig. 3b. Mean squared error was used as loss function for the patch-based regression head. The total loss for training our network was calculated by Eq. 1, where c specifies the ground truth grade, γ is a tuneable focusing parameter, α t the class imbalance weighting factor, p t is the model's estimated probability for the class with grade c = 1, and x,y are the arbitrary shapes. The network was trained with the Adam optimiser by using a maximal learning rate of 0.001 for 100 epochs until convergence was reached. Additionally, to minimise the number of anchors and therefore further optimise the architecture towards inference speed we only used the 32 × 32 feature map from the FPN. This was motivated by the fact that anchors of higher feature map sizes did not fit the small cell sizes and are limited in their total number. For comparison, we additionally tested Faster-RCNN 24 with a ResNet-50 backbone and SSD 25 with MobileNetV2 as provided by Huang et al. 35 . Both networks were trained with the Adam optimiser and a learning rate of 0.0001 for 100 epochs until convergence was reached. All networks were trained with random rotation, horizontal and vertical flips, but without intensity augmentations. This was appropriate since a shift in intensity could alter the cell grade.
Estimation based on image patch regression. Direct estimation of the hemosiderophages score by using an image patch-based regression approach is an alternative if the bounding box illustration is not required. Furthermore, an image patch-based regression approach could be used to find regions of interest efficiently even with standard computer vision approaches which we will discuss in the following two methods for a regression-based score estimation. While the first one used a support vector machine (SVM), the second was an adaptation of the RetinaNet architecture. The goal of the regression-based algorithm was to predict the grading score in a range the from zero to four on an image patch and to average the results for a total WSI. Support vector machine. In order to set a computationally inexpensive baseline for the task of estimating a hemosiderophages score, we trained a support vector machine with a Radial Basis Function (RBF) kernel and a convexity value of 0.1. These parameters were found by a grid search for the kernel and complexity parameter. As features we used the extracted histograms of a hundred patches per WSI with the sampling strategy described before.
Deep learning-based regression. To estimate the hemosiderophages score with a deep learning-based method we used the features extracted from RetinaNet and added two fully connected layers and a sigmoid activation for the regression head (see Fig. 3b). To compensate for the numerical instability of sigmoid activations close to zero and one and in order to enable a prediction score of up to grade four we scaled the sigmoid activation to a range from −0.5 and 4.5. The deep learning-based regression network was trained as a part of our RetinaNet-based object detection pipeline described in section object detection-based WSI score estimation.

Results
All experiments were run on a Linux workstation with a NVIDIA Quadro P5000 graphics card. The average calculation time for the object detection task was 101 seconds per WSI. The code for all experiments is available online and implemented in pytorch 40 27 and is commonly used to assess object detection performance. AP is the average detection precision under different recalls and mean Average Precision (mAP) is the average over all five grades.
Cell classification (CoCH). As stated above, we conducted an assessment of expert classification performance for comparison and to set a baseline. Comparing human experts and the deep learning classification pipeline, we found only offsets by one class by the deep learning system, whereas human expert disagreement (2020) 10:9795 | https://doi.org/10.1038/s41598-020-65958-2 www.nature.com/scientificreports www.nature.com/scientificreports/ was generally higher, especially for the higher grades 2, 3 and 4. In these categories disagreement in grade was significant for some cases (see Fig. 4). Concordance with the ground truth data was 85% for both automatic methods, whereas the human experts scored in a range of 69-86% (mean = 74, SD = 5) for the first round of labelling (V0) and 66-81% (mean = 73, SD = 4) for the second round of labelling (V1). This illustrates that we were able to reach human expert-level concordance with the cell-based regression and classification approach. The intra-observer variability ranges from 68 to 88% (mean = 79, SD = 6) with a mean Cohen's kappa score of 0.74. The inter-observer Fleiss' kappa score was 0.67 at the first round of labelling (V0) and 0. Object detection (MCWSA). Our object detection approach showed a mean average precision (mAP) of 0.66 (SD = 0.18, IoU = 0.5) over the three test set WSIs with a total of 3,518 patches and 15,427 cells. Table 2 shows the results per WSI over all tested networks with a maximal mAP of 0.66 reached by multiple approaches. The average error for cell-based grade score was 9 (SD = 24) and was calculated by taking the absolute difference of all ground truth cell grades and the predicted grades. For better understanding, a patch-wise analysed WSI is shown in Fig. 5.
The comparison of the three sampling strategies revealed a good overall convergence for the two stage cluster sampling strategy (mAP 0.66) and the quadtree sampling strategy (mAP 0.66), while completely random sampling showed very slow convergence to a lower mAP of 0.28.  www.nature.com/scientificreports www.nature.com/scientificreports/ Finally, this work has some limitations that need to be mentioned. All ground truth annotations were made by a single veterinary pathologist, data was collected at one laboratory and specimens were digitised with a single slide scanner. The data set comprises only seventeen WSIs, so our proposed approaches need to be validated on a larger, more diverse data set. Furthermore, we have taken no action to make an external colour calibration of the participants' screens which could positively influence the results of the participants but does not correspond to current clinical practice. In further work, we plan to analyse the effect of manual region selection by human experts and evaluate and reduce its impact on the proposed object detection pipeline. Furthermore, we are going to introduce an interactive annotation to increase the quality of the data set and effectively label new WSIs while analysing possible bias introduced by this. Also, it would be of high interest for us to analyse and evaluate our proposed methods and trained models on human pulmonary haemorrhage data sets. We believe a transfer to an application on human samples may be possible using only a small data set and transfer learning.

Data availability
The data sets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.