Risk stratification of malignant melanoma using neural networks
In order to improve the detection and classification of malignant melanoma,
this paper describes an image-based method that can achieve AUROC values of up
to 0.78 without additional clinical information. Furthermore, the importance of
the domain gap between two different image sources is considered, as it is
important to create usability independent of hardware components such as the
high-resolution scanner used. Since for the application of machine learning
methods, alterations of scanner-specific properties such as brightness,
contrast or sharpness can have strong (negative) effects on the quality of the
prediction methods, two ways to overcome this domain gap are discussed in this
paper.
Julian Burghoff, Leonhard Ackermann, Younes Salahdine, Veronika Bram, Katharina Wunderlich, Julius Balkenhol, Thomas Dirschka, Hanno Gottschalk
# Risk stratification of malignant melanoma using neural networks

###### Abstract

In order to improve the detection and classification of malignant melanoma, this paper describes an image-based method that can achieve AUROC values of up to 0.78 without additional clinical information. Furthermore, the importance of the domain gap between two different image sources is considered, as it is important to create usability independent of hardware components such as the high-resolution scanner used. Since for the application of machine learning methods, alterations of scanner-specific properties such as brightness, contrast or sharpness can have strong (negative) effects on the quality of the prediction methods, two ways to overcome this domain gap are discussed in this paper.

 neural nets, survival analysis, maligne melanom, domain gap

## 1 Introduction

Clinical and pathological staging of melanoma patients only relies on tumor size (Breslow thickness), ulceration and lymph node involvement [1]. However, the patient group with the thinnest melanomas in Tumor Stage T1 and with the most favorable prognosis resulted the most melanoma deaths in absolute numbers [2]. Additionally, patients from T3b onwards can be now offered potent adjuvant therapy [3, 4]. To identify patients with small tumors but high risk of relapse or to spare patients in advanced disease but with low mortality or recurrence risk, there is a current need for biomarkers and better prognostication [5]. The use of digitalized histological images like H&E scans, that help pathologists to better interpret the information provided by the tissue sample under the microscope, has been widely tried to use for improving diagnosis and prognosis in many tumors [6]. 3D high-resolution volumetric imaging of tissue architecture from large tissue and molecular structures at nanometer resolution are new techniques for improving early cancer detection, personalized risk assessment and potentially identifying the best treatment strategies [7]. Deep learning has been shown to read out additional information of these stains. These models have shown to deliver independent prognostic information[8, 9] and could also predict the results of molecular biomarkers [10].

In melanoma, CNN networks were shown to reach a concordance level above 80% for diagnosis compared to human pathologists[11], could outperform histopathologists in classification[12], and the result of the sentinel lymph node status [13].

Prognostication and predicting the risk of tumor recurrence could be shown for combining digital analysis with the detection of tumor-infiltrating lymphocytes by achieving a negative predicting value (NPV) of about 85% [14]. Another CNN approach for predicting disease specific survival in melanoma resulted in mixed results achieving area under the curve (AUROC) values of 90% and 88% but only NPVs of 95% and 65%, respectively, in two validation cohorts [15].

Such problems of the domain gap often arises in image recognition [16], i.e. a method that has been optimised on a data set from a certain source, but on data from a different source it provides unsatisfactory results. To address this problem, we first standardise the predictions for each dataset or data source, which already improves the accuracy. In a second approach, we add an additional regularisation term directly to the neural network so that the results are in the same range. This also improves the accuracy of the predictions.

The following parts of this paper are structured as follows: First, in section 2 we give a description of the data in terms of its clinical and technical structure. Section 3 provides an overview of methods that have been used in similar work and how our method is designed and built on these. Then, in section 4, we present and classify the results of our experiments. Finally, section 5 discusses the results and considers how future work might build on them.
## 2 Description of Data

In this paper we consider two main sets of data:

* _Dataset A_: This dataset contains 767 images of 176 patients of the American Joint Committee on Cancer (AJCC) set stages IA to IIID from four different locations in Bern (Switzerland), Bochum, Bonn and Kiel, Germany. All images have been recorded with the H&E colouring and were created by the same scanner Hamatasu NanoZoomer S210, NDP Version 2.4. In order to obtain a data set of the highest possible quality, the data were first manually checked by medical experts and a total of 23 data points were removed from the data set, e.g. due to broken slides. This clean-up took place before we split the data into training, validation and test data: The patients are assigned to train (104 patients), test (36 patients) and validation (36 patients) what results in 591 train-, 74 validation- and 102 test images. The split was chosen so that the distribution of high/low risk patients in each subset (training, validation, test) corresponds to the distribution of the full dataset A, see also table 1, under the constraint that multiple images of one patient remain in the same subset.
* _Dataset B_: The second dataset contains 242 images of 242 patients with AJCC stages IIA to IIC from the Central Malignant Melanoma Registry (CMMR) in Tubingen, Germany, where the H&E-colored images have partly different coloring than the images of _Dataset A_, probably caused by another scanner type Hamatasu Nanozoomer 2.0 HAT, NDP Version 2.5 and the scans are mostly disturbed by a marker pen on the slide. To see how our algorithm performs on images with a domain gap _dataset B_ is only used as a separate test dataset.

### Clinical Description

The images used represent hematoxylin-eosin (HE) stained melanoma sections. HE is a staining technique from histology that is used to better predict the disease prognosis of a melanoma patient, among other things, the mitosis rate can be determined (S3 Leitlmie Melanom).

Overall survival of dataset A (86%) is similar to dataset B (88%). However, the MSS time for dataset B with a median of 41 month was considerably lower compared to 70/ 73 months of dataset A. Relapse free survival differed as well with around 78% and 69% for dataset A and dataset B, respectively.

### Technical Description

Both datasets contain a total of 1009 images of different sizes (up to a resolution of \(158720\times 115456\) pixels) and various format ratios. Therefore, the images have to be pre-processed as neural nets on commercially available hardware are not yet able to handle images of this size at the time of writing. See section 3.2 for more information on the pre-processing task. The total size of the dataset is 570,3 GB in.ndpi file format.

For every patient \(i\), we have at least one image \(x_{i}\) and the information \(\delta_{i}=1\), indicating the death of the patient at time \(T_{i}\). If the patient survived the observation period of this study, we set \(\delta_{i}=0\) and \(T_{i}\) is the time the patient has been observed.

## 3 Materials and Methods

### Related Work

There exists abundant work on using deep neural networks for the interpretation of medical images. In this section, we focus on works with the scope of AI based prediction of survival for patients with a melanoma diagnosis.

In prior research on this topic [15] the authors use a two stage pipeline to predict risk on the basis of primary melanoma tumor images. Within this pipeline, they first use a segmentation to classify detect and crop tumor areas in the image. These small but detailed crops of \(500\times 5000\)px are then fed to a network of convolutional, recurrent and fully connected layers in order to predict the risk. [17] describes the approach to use a multivariable classifier that contains, besides clinical data, a score of a deep neural net, in order to predict the immunotherapy response of patients with advanced melanoma. Also here, clinical data is used for the model and a separation of segmentation and response classifier takes place.

The authors of [12] describe an experiment where a trained ResNet50 model outperforms 11 pathologists in classifying labeled histopathological images what shows that neural nets in general have high potential to improve correct melanoma diagnoses.

Similar to our approach, the authors of [18] used a VGG-based neural network architecture to detect cutaneous melanoma, although they use a binary classification for dead/alive patients instead of a survival analysis method. They evaluated their method a dataset provided by The Cancer Imaging Archive of 53 patients with a given survival status.

In [19] the authors on the one hand locate molecular biomarkers in immunohistochemistry images using convolutional neural networks which can help enabling new cancer screenings. On the other hand they also classify the found biomarkers along their type.

However, there are not only machine learning approaches to melanoma classification based on H&E scans, but also based on topics such as dermoscopic image data [20] where the authors use a quiet small 5-layers CNN to classify the images to the corresponding tumor stage with applying the Adam optimizer on the Similarity Measure for Text Processing as loss function. This work is also based on [21] where the authors predict the melanomas thickness using a pretrained VGG-19 model on \(400\times 400\)px preprocessed dermoscopic images.
In [22] the authors also use deep convolutional neural networks in combination with survival analysis in order to find good predictions on pathological images - here in context of lung cancer. They annotated the regions of interest of the images with help of pathologists and sampled small random high resolution crops of these regions to use in the networks whereas we used down-sampled low resolution images in our approach.

In contrast to this, in our work we extract information directly from the original image data using a VGG16-like neural network, so we don't use an additional pre-segmentation, which might be error prone by itself. In this way we also avoid labeling of the regions of interest by humans annotators, which is a time consuming process and requires highly trained annotators.

### Preprocessing the Data

When pre-processed, the images are reduced from their original format by a factor of 64 in dataset A and 128 in dataset B in each dimension. The resulting image is centred in a 2500 by 2000 pixel frame which is filled with white color outside the image.

Images of patients with multiple images are seen as independent information in the training data.

### Methods

#### 3.3.1 Cox's Proportional Hazards Model

Survival analysis is about predicting the probability of absence of an event (e.g. the death of a patient) until time \(t\) using the parameters \(\beta\) of a suitable model [23]. One of the most widely used methods is Cox's Proportional Hazards Model, which predicts the hazard function \(h(t|x)\) on the basis of an input vector \(x\), see e.g., [15, 22] or [24] :

\[h(t|x)=h_{0}(t)\cdot\exp(\beta^{T}x)\Leftrightarrow\log\frac{h(t|x)}{h_{0}(t) }=\beta^{T}x \tag{1}\]

The baseline hazard function \(h_{0}(t)\) indicates how large the hazard rate would be without the influence of other parameters (like \(\beta\)) and is therefore only dependent on the time \(t\), which means that for our case it is eliminated from the calculation of the loss function and therefore does not need to be calculated [24]. To estimate the parameters \(\beta\) of the linear model the negative log partial likelihood function can be minimized:

\[l(\beta)=-\sum_{i=1}^{n}\delta_{i}(\beta^{T}x_{i}-\log\sum_{j\in R_{i}}\exp( \beta^{T}x_{j})) \tag{2}\]

where \(\delta\) is the delta-function which determines whether the data is censored or not, \(n\) the total number of data points and \(R_{i}=\{j|y_{j}\geq y_{i}\}\) is the risk set which describes the data-subset of patients which do not had an event before timestamp \(y_{i}\) of the \(i\)-th event.

An adaption of linear models to non-linear models like (deep) neural nets was introduced by [25] in 1995 and applications can be found for example in [22] or [24]. The risk function \(\beta^{T}x\) is replaced by the output of the neural net \(\hat{h}_{\theta}(x)\) and we achieve the following loss function for each mini batch \(B\):

\[\mathcal{L}(\theta)\coloneqq-\sum_{i\in B}\delta_{i}\left(\hat{h}_{\theta}(x_ {i})-\log\sum_{j\in R_{i}}e^{\hat{h}_{\theta}(x_{j})}\right) \tag{3}\]

#### 3.3.2 Neural Net Architecture and Training

Convolutional neural networks (CNN) represent the state of the art in image recognition. As neural network model, we use a modified version of a VGG16 net [26] in our experiments. Figure 1 gives an overview of the network architecture. Each convolutional layer has kernel sizes of \(3\times 3\) and each pooling layer is a maximum pooling with pooling size \(2\times 2\).

As mentioned in section 3.1, we use the loss function supplied in [24] which implements a variant of the Cox's Proportional Hazards Model. In order to achieve a higher generalisability and reducing the domain gap between dataset A and dataset B, we expanded our loss function in 3 by a regularization term:

\[\mathcal{L}(\theta)\coloneqq-\sum_{i\in B}\delta_{i}(\hat{h}_{ \theta}(x_{i})-\log\sum_{j\in R_{i}}e^{\hat{h}_{\theta}(x_{j})})\\ +\lambda(\frac{1}{|B|}(\sum_{j\in B}\hat{h}_{\theta}(x_{j}))^{2} \tag{4}\]

where \(\lambda\) is the chosen regularization strength. See sec 4.1 for further information why we use this additional regularization term.

\begin{table}
\begin{tabular}{|c|c c c|c|} \hline Characteristic & \multicolumn{3}{c|}{Dataset A} & \multicolumn{1}{c|}{Dataset B} \\  & Training & Validation & Test & Test \\ \hline N(Scans) & 313 & 36 & 38 & 242 \\ N(Patients) & 104 & 36 & 36 & 242 \\ \hline Alive / Censored & 89 (86\%) & 31 (86\%) & 31 (86\%) & 212 (88\%) \\ Dead / Event & 15 (14\%) & 5 (14\%) & 5 (14\%) & 30 (12\%) \\ MSS time (months) Mean & 70.12 & 80.03 & 78.37 & 50 \\ MSS time (months) Median & 70 & 73 & 73 & 41 \\ \hline Relapse Free Survival Recurrence & 21 (20\%) & 8 (22\%) & 7 (19\%) & 75 (31\%) \\ Relapse Free Survival Non-recurrence & 83 (80\%) & 28 (78\%) & 29 (81\%) & 167 (69\%) \\ RFS time Mean & 67 & 76.67 & 75 & 42 \\ RFS time Median & 70 & 68 & 73 & 30 \\ \hline \end{tabular}
\end{table} TABLE 1: Important features of our datasets.
In training, we employ the Adam optimizer over \(50\) epochs with a learning rate of \(\alpha=0.001\) and a batch size of \(5\).

For the implementation we used Openslide (version 1.1.2) [27] for importing the images, the GitHub repository of Sebp[24] and Tensorflow (version 2.6.1) [28] along with Keras (version 2.6.0) [29].

For the calculations, we used a workstation with a Dual Intel Xeon Gold 6248R 3.0GHz and three Nvidia Quadro RTX 8000 graphic units with 48GB VRAM each, whereas the different GPUs are only used for different trainings.

#### 3.3.3 Concordance Index

For validation purposes of our methods we use Harrel's concordance index (CI) [30]. The CI is defined as the ratio between correctly ordered pairs and all possible rankable pairs [31]:

\[CI=\frac{\#\text{concordant pairs}}{\#\text{comparable pairs}} \tag{5}\]

A pair of observations \(i,j\) with its survival times fulfill \(T_{i}>T_{j}\), is concordant if \(\hat{h}_{\theta}(x_{j})>\hat{h}_{\theta}(x_{i})\). Also a pair \(i,j\) is not comparable if the smaller survival time is censored (i.e. \(T_{i}>T_{j}\wedge\Delta_{j}=0\)). Otherwise this pair is comparable. Thus,

\[CI=\frac{\sum_{i,j}\mathds{1}(T_{i}>T_{j})\cdot\mathds{1}(\hat{h}_{\theta}( x_{j})>\hat{h}_{\theta}(x_{i}))\cdot\Delta_{j}}{\sum_{i,j}\mathds{1}(T_{i}>T_ {j})\cdot\Delta_{j}} \tag{6}\]

The CI estimates the probability of concordance \(P(\hat{h}_{\theta}(x_{j})>\hat{h}_{\theta}(x_{i})|T_{i}>T_{j})\) for two independent observations/predictions. It can also be interpreted as a measure of the area under a time-dependent receiver operator curve [32][31],[33]. A value of \(CI=1\) means that all observations are correctly sequenced, \(CI=0.5\) means that the method applied is no better than guessing.

#### 3.3.4 Area Under the Receiver Operating Characteristic

In the evaluation of the results our methods generate an important measurement value is the area under the receiver operating characteristic (AUROC)[34]. For a data point with \(a=P(F_{p})\) the probability of the data point being a false positive prediction and \(b-1=P(T_{p})\) the negative probability of the data point being a true positive prediction in a dataset \(D\), [35] defines the AUROC as follows:

\[\mathrm{AUROC}=\sum_{i\in D}\{(1-b_{i}\cdot\Delta a)+\frac{1}{2}[\Delta(1-b) \cdot\Delta a]\}\]

with

\[\Delta(1-b)=(1-b_{i})-(1-b_{i-1})\]

and

\[\Delta a=a_{i}-a_{i-1}\]

## 4 Results

An overview of all the concordance indices and AUROC values of our experiments is provided in table 2.

### Results of model without regularization

In our first set of experiments, we use the loss function supplied in equation (3). Although we achieve an AUROC value of \(61.5\%\) on the test subset of dataset A and also an AUROC value of \(63.5\%\) on the separated dataset B. We observe that the domains in which our predictions \(\hat{h}_{\theta}(x_{i})\) lay differ significantly from dataset A to dataset B.

This results in a bad overall AUROC value of \(54.4\%\) if one mixes these predictions naivly together as it would be a complete dataset from only one datasource before evaluating the whole set. One potential reason is that the net learns features that are relevant for the task, but it also is sensitive to further properties of the image like the pixel's brightness (or in general the pixel's color distribuition). While such image features do not encode meaningful medical information, they can still disturb the outcome of the network, especially if the

Figure 1: Design of used neural network. Orange are convolutional layers, red are pooling layers and purple indices dense layers.
hazards predicted on the second data set are in a different numerical range as compared with the original training data. This in particular happens by deviation from the neutral direction \(\hat{h}_{\theta}(x_{j})\to\hat{h}_{\theta}(x_{j})+z\) which merely leads to a redefinition of the baseline hazard function \(h_{0}(t)\) but is not sensed in the Cox loss function. So one can imagine the predicitons \(\hat{h}_{\theta}(x_{i})\) shifted away from the total diagonal like schematically shown in figure 2.

In our case, dataset B has a slightly other color scheme (mainly because of the use of another scanner - see section 2) and so our predictions evaluating the model trained on dataset A on this dataset leads to a significant drop in performance.

A first approach to reduce the sensitivity towards different hazards is to re-center the predicted hazards. We thus standardise (we subtract the mean and divide the result by the standard deviation) the predictions of each dataset and merge them afterwards. We display the resulting improvement in the rightbest column of each approach in table 2.

The downside of this approach is that you need to have already a certain set of of images for each source of images/scanner that is used to predict survival probabilities.

### Results of model with regularization

To interdict shifting of hazard functions altogether, we add a L2-regularization-term, see equation (4), to shift all predictions \(\hat{h}_{\theta}(x_{i})\) in the same domain range. This does not only lead to a higher AUROC value in evaluation on the mixed data set (\(60.1\%\) instead of \(54.4\%\)), but we can improve even more if combining the regularization with the above-mentioned normalization process (\(64.2\%\) instead of \(61.4\%\))

Besides, we also achieved a significantly higher ROC value on the testdata of dataset A throughout this regularization technique. Unfortunately the overall generalizability on the isolated dataset B suffered (\(57.8\%\) instead of \(63.5\%\)) from that approach.

## 5 Discussion and Outlook

We have shown that it is possible to make survival predictions based on simplified image information using Cox's Proportional Hazard Models on neural networks and that our domain adaptation techniques succeed in merging the predictions into one range. Compared to [18] where the authors achieved an AUROC value of 76.9% our approach with regularization slightly outperforms it with an AUROC value of 79.5%.

In future work, we will have a more detailed evaluation and we have a closer look on levels of significance and correlations to clinical variables like the Breslow depth or genetic information and demographic variables such as age or gender to evaluate and improve the clinical utility of our predictions.
