# COMBINING HIGH-LEVEL FEATURES OF RAW AUDIO WAVES AND MEL-SPECTROGRAMS FOR
AUDIO TAGGING

###### Abstract

In this paper, we describe our contribution to Task 2 of the DCASE 2018 Audio
Challenge [1]. While it has become ubiquitous to utilize an ensemble of
machine learning methods for classification tasks to obtain better predictive
performance, the majority of ensemble methods combine predictions rather than
learned features. We propose a single-model method that combines learned high-
level features computed from log-scaled mel-spectrograms and raw audio data.
These features are learned separately by two Convolutional Neural Networks,
one for each input type, and then combined by densely connected layers within
a single network. This relatively simple approach along with data augmentation
ranks among the best two percent in the Freesound General-Purpose Audio
Tagging Challenge on Kaggle.

Index Terms—  audio-tagging, convolutional neural network, raw audio, mel-
spectrogram

##  1 INTRODUCTION

For humans, it seems to be effortless to associate sounds with events or
categories that describe the perceived sound best. However, the complex
structure and the large amount of information transmitted through sound makes
it particularly difficult to extract that information automatically.

Recognizing a wide variety of sounds has many applications in our today’s
life. These include surveillance [2, 3, 4], acoustic monitoring [5], and
automatic description of multimedia [6]. Due to the diversity of sounds
belonging to the same category, a reliable recognition of manifold sound
categories is still under ongoing research.

Carefully hand-crafted features such as Mel Frequency Cepstral Coefficients
(MFCCs) were the dominant features used for speech recognition [7, 8] and
music information retrieval [9], but the trend is now shifting toward deep
learning [10]. The approach of manual feature engineering has drawbacks
compared to deep learning based methods because it requires considerable
effort and expertise to manually create features for a specific purpose. In
particular, most of the engineered features, such as MFCCs and spectral
centroids [11], are non-task specific, whereas all deep learning approaches
are task-specific due to its formulation as a minimization process on task-
specific training examples. Since Convolutional Neural Networks (CNNs) have
shown a remarkable progress in visual recognition tasks [12] over the last
years, it has become common to use CNNs for feature extraction and
classification in the audio domain [13, 14]. Several CNN architectures, such
as AlexNet [15], VGG [16], ResNet [17], and Inception-v3 [18], have been
proposed for image classification, which are also well suited to the task of
audio tagging.

The goal of Task 2 of the DCASE 2018 Challenge [1] was to predict the category
of an audio clip belonging to one out of 41 heterogeneous classes, such as
“Acoustic guitar”, “Bark”, “Bus” and “Telephone”, drawn from the AudioSet
Ontology [19]. Training and testing data contain a diverse set of user-
generated audio clips from Freesound (https://freesound.org) [1].

In this paper, we mainly focus on building an audio-tagging system that uses
both the raw audio data and the corresponding mel-spectrogram rather than
ensembling [20] or stacking [21] multiple classifiers.

##  2 METHOD

Our audio-tagging system comprises two separately trained Convolutional Neural
Networks on raw audio and mel-spectrogram, respectively. The learned high-
level features are then combined by a densely connected neural network to form
the system. In the following, we describe each model in detail.

###  2.1 CNN on Raw Audio (cnn-audio)

For the cnn-audio model, we use an architecture similar to common
architectures for image classification, like VGG16 [16] or AlexNet [15], but
with one-dimensional convolutions and one-dimensional max pooling. As
described in Table 1, we use four blocks, each consisting of two convolutional
and one max-pooling layer. The number of filters is increased in each
consecutive block, while the kernel size is decreased. The pool-size of the
max-pooling layers is chosen to quickly reduce the large time dimension. After
each block and after the dense layer, we apply batch normalization [22], as
experiments have shown that it reduces the training time and increases the
model accuracy. To introduce nonlinearities into our network, we apply a ReLU
activation function [23] after each convolutional layer and dense layer.

Layer | 1sec shape | 2sec shape | 3sec shape
---|---|---|---
input | (44100, 1) | (88200, 1) | (132300, 1)
conv1d, 11, 32 | (44100, 32) | (88200, 32) | (132300, 32)
conv1d, 11, 32 | (44100, 32) | (88200, 32) | (132300, 32)
max-pool1d, 8/16 | (5512, 32) | (5512, 32) | (8268, 32)
conv1d, 9, 64 | (5512, 64) | (5512, 64) | (8268, 64)
conv1d, 9, 64 | (5512, 64) | (5512, 64) | (8268, 64)
max-pool1d, 16 | (344, 64) | (344, 64) | (516, 64)
conv1d, 7, 128 | (344, 128) | (344, 128) | (516, 128)
conv1d, 7, 128 | (344, 128) | (344, 128) | (516, 128)
max-pool1d, 16 | (21, 128) | (21, 128) | (32, 128)
conv1d, 5, 256 | (21, 256) | (21, 256) | (32, 256)
conv1d, 5, 256 | (21, 256) | (21, 256) | (32, 256)
max-pool1d, 16 | (1, 256) | (1, 256) | (2, 256)
dense, 512 | (512) | (512) | (512)
softmax, 41 | (41) | (41) | (41)
Table 1: The architecture of the cnn-audio model. Note that for the one-second
model, the first max-pooling layer uses a pool-size of 8, while for other
models a pool-size of 16 is used.

###  2.2 CNN on Mel-Spectrogram (cnn-spec)

The cnn-spec model is a two-dimensional convolutional neural network taking
mel-spectrograms as input. The architecture is again similar to common image
classification architectures and is described in detail in Table 2. As with
the one-dimensional model, we apply batch normalization after each block and
the dense layer and use the ReLU activation function after convolutional and
dense layers.

The mel-spectrogram is extracted using librosa [24] with the original sampling
frequency of $44.1$ kHz, 2048 FFT points, 128 mel-bins, and a hop-length of
256. The amplitude of the mel-spectrogram is scaled logarithmically, and the
scaled mel-spectrogram is resized in time dimension to fit the model input
size.

Layer | 1sec shape | 2sec shape | 3sec shape
---|---|---|---
input | (128, 170, 1) | (128, 300, 1) | (128, 400, 1)
conv2d, 4$\times$4, 64 | (128, 170, 64) | (128, 300, 64) | (128, 400, 64)
conv2d, 4$\times$4, 64 | (128, 170, 64) | (128, 300, 64) | (128, 400, 64)
max-pool2d, 1$\times$1/2$\times$2 | (128, 170, 64) | (64, 150, 64) | (64, 200, 64)
conv2d, 4$\times$4, 64 | (128, 170, 64) | (64, 150, 64) | (64, 200, 64)
max-pool2d, 2$\times$2 | (64, 85, 64) | (32, 75, 64) | (32, 100, 64)
conv2d, 3$\times$3, 128 | (64, 85, 128) | (32, 75, 128) | (32, 100, 128)
max-pool2d, 2$\times$4 | (32, 21, 128) | (16, 18, 128) | (16, 25, 128)
conv2d, 3$\times$3, 128 | (32, 21, 128) | (16, 18, 128) | (16, 25, 128)
max-pool2d, 2$\times$2 | (16, 10, 128) | (8, 9, 128) | (8, 12, 128)
conv2d, 3$\times$3, 256 | (16, 10, 256) | (8, 9, 256) | (8, 12, 256)
max-pool2d, 2$\times$2 | (8, 5, 256) | (4, 4, 256) | (4, 6, 256)
conv2d, 3$\times$3, 256 | (8, 5, 256) | (4, 4, 256) | (4, 6, 256)
max-pool2d, 2$\times$2 | (4, 2, 256) | (2, 2, 256) | (2, 3, 256)
dense, 256 | (256) | (256) | (256)
softmax, 41 | (41) | (41) | (41)
Table 2: The architecture of the cnn-spec model. Note that the first max-
pooling layer does not exist in the case of the one-second model.

###  2.3 Joining CNNs (cnn-comb)

We remove the softmax and dense layer of both the trained cnn-audio and cnn-
spec model and then concatenate the output features of the previous layer of
both models so as to join them. The concatenated features are then connected
to a densely connected neural network with four hidden layers. The hidden
dense layers have 512, 256, 256, and 128 neurons, respectively. The complete
model is illustrated in Figure 1.

We train cnn-audio and cnn-spec from scratch. Afterward, the weights of these
models are transferred to the cnn-comb model and only the newly added dense
layers are trained. The splitting of the training of cnn-comb into three
steps, facilitates the procedure.

Figure 1: Illustrated architecture of the complete model.

###  2.4 Data Augmentation

To prevent our model from overfitting, we make use of extensive data
augmentation during training (time shifting, cropping, padding, and blending
clips of same and different categories). Each of these augmentation techniques
is applied to the raw audio wave and the mel-spectrogram. In the remaining
section, we explain the augmentation methods based on the raw audio wave.

First, we apply a uniformly random time shift to the audio clip. To ensure
that the audio clips fit the size of the model input, crops are taken from too
long audio files and too short audio files are padded. For audio clips that
are longer than the model input size, we use a crop with the size of the model
input taken from a random position. If an audio sample fits multiple times
($n$ times) in the input size, it is replicated such that it appears
$k\in\\{1,\ldots,n\\}$ times with a probability of $1/n$, and the remaining
space before, after, and in between replications is filled with zeros.

Additionally, we enhance the robustness of our model by blending multiple
audio clips of same or different categories. This method is referred to as
mixup [25]. We blend them by assigning a random weight to each sample (weights
sum up to one) and taking the weighted sum. If the blended samples are of the
same class, the model should still predict the common class for the newly
generated training sample. In the other case, if the blended samples are of
distinct classes, the model is trained to predict the weight of each included
class.

While the spectrogram is computed in advance, the computationally inexpensive
data augmentation techniques can be computed on-the-fly during training. This
saves disk-space and guarantees a large amount of diverse training data.

###  2.5 Implementation Details

We implemented the described method using Keras [26] in Python.

To monitor overfitting and the model performance during training, we exclude a
part of the training data as validation data. To still make use of all
training data, we train five models on stratified folds of the training data
such that each training example is used once for validation and four times for
training. For the final prediction, we accumulate the predictions of all five
models using the geometric mean.

All models are trained using the Adam optimizer [27] with a fixed learning
rate of $0.001$ and a mini-batch size of $32$ for a maximum of $300$ epochs,
but stopping earlier if the validation loss hasn’t improved for $35$ epochs.
We use the categorical cross-entropy loss function and weight the loss
according to the distribution of training examples per class, thereby ensuring
that the models pay more attention to samples from an under-represented class.

Because many clips contain silence, we cut off silent parts at the beginning
and at the end of an audio clip that do not exceed a volume of $40$ decibel.

When predicting the test data, we have to take the varying length of audio
files into account. It is not sufficient to only predict on one crop of too
long tracks because important features might not be present in the selected
crop. Therefore, we run the inference on many crops of the audio file with a
step size of $5120$ frames, which is approximately $0.12$ seconds. For too
short audio tracks, the model might be able to better recognize class-specific
features in certain parts of the input. Therefore, we generate multiple inputs
by padding the audio file with zeros such that the real audio appears at
different positions in the input. Again, we use a step size of $5120$ frames.
Multiple predictions for one audio file are combined by means of the geometric
mean.

##  3 EVALUATION

###  3.1 Dataset

We evaluate our method on the dataset provided for Task 2 of the DCASE 2018
Challenge [1], which comprises $9473$ training and $1600$ test samples. The
test data set has been manually verified, whereas the training data features
labels of different reliability. Each mono audio file has a bit-depth of $16$,
a sampling rate of $44.1$ kHz, and is associated with one out of $41$ classes
of the AudioSet Ontology [19].

The class distribution of both training and test set is not balanced and
ranges from $94$ to $300$ and from $25$ to $110$ samples per class,
respectively. The duration of the shortest audio file is $300$ ms and $30$
seconds for the longest clip, while the average length is $6.8$ seconds for
the train set and $5.2$ seconds for the test set.

###  3.2 Metric

The Challenge uses mean Average Precision at three (mAP@3) for evaluating test
results, which allows up to three predictions per audio clip. Full credit is
given if the first prediction matches the label of the clip, while less credit
is given if one of the other predictions is correct. The evaluation metric is
defined as

$\text{mAP@3}=\frac{1}{U}\sum_{i=1}^{U}\sum_{j=1}^{\min(3,n_{i})}\frac{\llbracket
y_{ij}=\hat{y}_{i}\rrbracket}{j},$

where $U$ is the total number of scored audio files, $y_{ij}$ is the predicted
label for file $i$ at position $j$, $\hat{y}_{i}$ is the ground-truth label
for file $i$, $n_{i}$ is the total number of predicted labels for file $i$,
$\llbracket\texttt{True}\rrbracket=1,\llbracket\texttt{False}\rrbracket=0$. No
label may be predicted multiple times for one audio file.

###  3.3 Results

We have trained all models on inputs of one, two, and three seconds, as
described in Section 2.5, and evaluated the model performance on the test set
(see Table 3). We have observed that for each crop length, the combined model
performs significantly better compared to models with a single input. Combined
models with an input size of two and three seconds, perform best and rank in
the upper two percent on the Private Leaderboard on Kaggle.

Model |  Crop length |  Public score (mAP@3) |  Private score (mAP@3) |  Total (mAP@3)
---|---|---|---|---
cnn-audio |  1sec |  0.920 |  0.888 |  0.894
2sec |  0.921 |  0.884 |  0.891
3sec |  0.935 |  0.889 |  0.898
cnn-spec |  1sec |  0.930 |  0.923 |  0.924
2sec |  0.950 |  0.928 |  0.932
3sec |  0.935 |  0.930 |  0.931
cnn-comb |  1sec |  0.955 |  0.939 |  0.942
2sec |  0.966 |  0.944 |  0.948
3sec |  0.956 |  0.944 |  0.946
Table 3: Evaluation results of the individual models on the public ($301$
samples), private ($1299$ samples), and full test set.

Additionally, we determined the per-category mAP@3 on the complete test set
showing that some classes are more challenging to predict than others (see
Table 4). Our model primarily struggles with the classes “Squeak”, “Telephone”
and “Fireworks”, but it still beats the baseline system [1] in every per-
category mAP@3 score.

Name | samples | time | mAP@3 |  Name | samples | time | mAP@3 |  Name | samples | time | mAP@3
---|---|---|---|---|---|---|---|---|---|---|---
Acoustic_guitar | 300 | 52.2 | 0.893 |  Electric_piano | 150 | 25.5 | 1.000 |  Microwave_oven | 146 | 25.1 | 0.966
Applause | 300 | 58.2 | 1.000 |  Fart | 300 | 18.6 | 0.944 |  Oboe | 299 | 15.3 | 0.976
Bark | 239 | 44.6 | 0.982 |  Finger_snapping | 117 | 5.9 | 1.000 |  Saxophone | 300 | 33.7 | 0.942
Bass_drum | 300 | 12.8 | 1.000 |  Fireworks | 300 | 48.2 | 0.786 |  Scissors | 95 | 15.7 | 0.927
Burping_or_eructation | 210 | 11.7 | 1.000 |  Flute | 300 | 46.2 | 1.000 |  Shatter | 300 | 26.1 | 0.960
Bus | 109 | 28.4 | 0.953 |  Glockenspiel | 94 | 8.4 | 0.856 |  Snare_drum | 300 | 17.9 | 0.912
Cello | 300 | 37.3 | 0.951 |  Gong | 292 | 41.8 | 0.968 |  Squeak | 300 | 38.2 | 0.603
Chime | 115 | 23.8 | 0.891 |  Gunshot_or_gunfire | 147 | 11.1 | 0.950 |  Tambourine | 221 | 10.1 | 0.975
Clarinet | 300 | 34.7 | 0.991 |  Harmonica | 165 | 18.6 | 0.970 |  Tearing | 300 | 38.7 | 0.981
Computer_keyboard | 119 | 23.0 | 1.000 |  Hi-hat | 300 | 18.6 | 0.957 |  Telephone | 120 | 16.2 | 0.788
Cough | 243 | 22.4 | 1.000 |  Keys_jangling | 139 | 18.8 | 0.929 |  Trumpet | 300 | 28.3 | 0.959
Cowbell | 191 | 10.9 | 1.000 |  Knock | 279 | 19.6 | 0.957 |  Violin_or_fiddle | 300 | 26.6 | 0.986
Double_bass | 300 | 16.9 | 0.946 |  Laughter | 300 | 36.3 | 0.974 |  Writing | 270 | 48.3 | 0.948
Drawer_open_or_close | 158 | 18.0 | 0.925 |  Meow | 155 | 18.7 | 1.000 |  |  |  |
Table 4: Per category mAP@3 score of the cnn-comb 2sec model on the full test
set and the number of samples along with time in minutes of the respective
class in the train set.

To verify that the performance gain of the combined model results from
combining the extracted high-level features from both models, we compare cnn-
audio and cnn-spec to the cnn-comb model, where one of its inputs is set to
zero (See Figure 2).

For categories in which the single cnn-audio model outperforms the single cnn-
spec model, the combined model performs better if the audio input is present
and not set to zero. Otherwise, if cnn-spec outperforms cnn-audio, cnn-comb
has a higher score if the mel-spectrogram is present and not set to zero.
Setting either the mel-spectrogram or the audio wave to zero forces the cnn-
comb model to make predictions based on a single input.

cnn-comb performs equally or better than cnn-comb with one single input set to
zero because it has learned to utilize meaningful high-level features of both
inputs jointly, which are not given by zeroed inputs. For the same reason,
ccn-comb with one zeroed input usually performs worse than the corresponding
model with a single input.

We conclude that cnn-comb makes use of the learned high-level features from
both cnn-audio and cnn-spec, but it focuses more on features belonging to the
superior single model.

Figure 2: Comparison of per-category scores of single input models, combined
models with one input alternately set to zero, and the combined model with
both inputs. The mAP@3 score is reported on a single fold for each model.

##  4 CONCLUSION

In this paper, we have proposed a method for audio-tagging that extends
current Convolutional Neural Network approaches that only make use of a
frequency representation by adding a second input that incorporates the raw
audio wave. Adding the additional input, has improved the mAP@3 score
significantly. We have demonstrated the capabilities of our model by competing
in the Freesound General-Purpose Audio Tagging Challenge on Kaggle and ranking
in the top two percent of all participants.

##  5 ACKNOWLEDGMENT

We thank Christian Borgelt and Christoph Doell for motivating us to take part
in the Kaggle competition and Christoph Doell for his valuable comments on the
manuscript.

## References

  * [1] E. Fonseca, M. Plakal, F. Font, D. P. Ellis, X. Favory, J. Pons, and X. Serra, “General-purpose tagging of freesound audio with audioset labels: Task description, dataset, and baseline,” _arXiv preprint arXiv:1807.09902_ , 2018\.
  * [2] A. Harma, M. F. McKinney, and J. Skowronek, “Automatic surveillance of the acoustic activity in our living environment,” in _Multimedia and Expo, 2005\. ICME 2005. IEEE International Conference on_. IEEE, 2005, pp. 4–pp.
  * [3] M. Crocco, M. Cristani, A. Trucco, and V. Murino, “Audio surveillance: a systematic review,” _ACM Computing Surveys (CSUR)_ , vol. 48, no. 4, p. 52, 2016.
  * [4] G. Valenzise, L. Gerosa, M. Tagliasacchi, F. Antonacci, and A. Sarti, “Scream and gunshot detection and localization for audio-surveillance systems,” in _Advanced Video and Signal Based Surveillance, 2007. AVSS 2007. IEEE Conference on_. IEEE, 2007, pp. 21–26.
  * [5] S. Goetze, J. Schroder, S. Gerlach, D. Hollosi, J.-E. Appell, and F. Wallhoff, “Acoustic monitoring and localization for social care,” _Journal of Computing Science and Engineering_ , vol. 6, no. 1, pp. 40–50, 2012.
  * [6] Q. Jin and J. Liang, “Video description generation using audio and visual cues,” in _Proceedings of the 2016 ACM on International Conference on Multimedia Retrieval_. ACM, 2016, pp. 239–242.
  * [7] C. Ittichaichareon, S. Suksri, and T. Yingthawornsuk, “Speech recognition using MFCC,” in _International Conference on Computer Graphics, Simulation and Modeling (ICGSM’2012) July_ , 2012, pp. 28–29.
  * [8] B. Logan _et al._ , “Mel frequency cepstral coefficients for music modeling.” in _ISMIR_ , vol. 270, 2000, pp. 1–11.
  * [9] G. Tzanetakis and P. Cook, “Musical genre classification of audio signals,” _IEEE Transactions on speech and audio processing_ , vol. 10, no. 5, pp. 293–302, 2002.
  * [10] O. Abdel-Hamid, L. Deng, and D. Yu, “Exploring convolutional neural network structures and optimization techniques for speech recognition.” in _Interspeech_ , vol. 2013, 2013, pp. 1173–5.
  * [11] A. J. Eronen, V. T. Peltonen, J. T. Tuomi, A. P. Klapuri, S. Fagerlund, T. Sorsa, G. Lorho, and J. Huopaniemi, “Audio-based context recognition,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 14, no. 1, pp. 321–329, 2006.
  * [12] O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein, _et al._ , “ImageNet large scale visual recognition challenge,” _International Journal of Computer Vision_ , vol. 115, no. 3, pp. 211–252, 2015.
  * [13] O. Abdel-Hamid, A.-r. Mohamed, H. Jiang, L. Deng, G. Penn, and D. Yu, “Convolutional neural networks for speech recognition,” _IEEE/ACM Transactions on audio, speech, and language processing_ , vol. 22, no. 10, pp. 1533–1545, 2014.
  * [14] K. J. Piczak, “Environmental sound classification with convolutional neural networks,” in _Machine Learning for Signal Processing (MLSP), 2015 IEEE 25th International Workshop on_. IEEE, 2015, pp. 1–6.
  * [15] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “ImageNet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [16] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _arXiv preprint arXiv:1409.1556_ , 2014.
  * [17] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [18] C. Szegedy, V. Vanhoucke, S. Ioffe, J. Shlens, and Z. Wojna, “Rethinking the inception architecture for computer vision,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 2818–2826.
  * [19] J. F. Gemmeke, D. P. Ellis, D. Freedman, A. Jansen, W. Lawrence, R. C. Moore, M. Plakal, and M. Ritter, “Audio Set: An ontology and human-labeled dataset for audio events,” in _Acoustics, Speech and Signal Processing (ICASSP), 2017 IEEE International Conference on_. IEEE, 2017, pp. 776–780.
  * [20] L. Breiman, “Bagging predictors,” _Machine learning_ , vol. 24, no. 2, pp. 123–140, 1996.
  * [21] S. Džeroski and B. Ženko, “Is combining classifiers with stacking better than selecting the best one?” _Machine learning_ , vol. 54, no. 3, pp. 255–273, 2004.
  * [22] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deep network training by reducing internal covariate shift,” _arXiv preprint arXiv:1502.03167_ , 2015.
  * [23] V. Nair and G. E. Hinton, “Rectified linear units improve restricted boltzmann machines,” in _Proceedings of the 27th international conference on machine learning (ICML-10)_ , 2010, pp. 807–814.
  * [24] B. McFee, M. McVicar, S. Balke, C. Thomé, C. Raffel, O. Nieto, E. Battenberg, D. Ellis, R. Yamamoto, J. Moore, R. Bittner, K. Choi, F.-R. Stöter, S. Kumar, S. Waloschek, Seth, R. Naktinis, D. Repetto, C. F. Hawthorne, C. Carr, hojinlee, W. Pimenta, P. Viktorin, P. Brossier, J. F. Santos, JackieWu, Erik, and A. Holovaty, “librosa/librosa: 0.6.0,” Feb. 2018. [Online]. Available: https://doi.org/10.5281/zenodo.1174893
  * [25] H. Zhang, M. Cisse, Y. N. Dauphin, and D. Lopez-Paz, “mixup: Beyond empirical risk minimization,” in _International Conference on Learning Representations_ , 2018.
  * [26] F. Chollet _et al._ , “Keras,” https://keras.io, 2015.
  * [27] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.


heightadjust=object capbtabboxtable[][]

# Continual atlas-based segmentation of prostate MRI

Amin Ranem1
Camila González2
Daniel Pinto dos Santos3
Andreas M. Bucher4
Ahmed E. Othman5
Anirban Mukhopadhyay1
1 Technical University of Darmstadt, 2 Stanford University, 3 University of
Cologne,
4 University of Frankfurt, 5 University Medical Center Mainz
<EMAIL_ADDRESS>
###### Abstract

Continual learning (CL) methods designed for natural image classification
often fail to reach basic quality standards for medical image segmentation.
Atlas-based segmentation, a well-established approach in medical imaging,
incorporates domain knowledge on the region of interest, leading to
semantically coherent predictions. This is especially promising for CL, as it
allows us to leverage structural information and strike an optimal balance
between model rigidity and plasticity over time. When combined with privacy-
preserving prototypes, this process offers the advantages of rehearsal-based
CL without compromising patient privacy. We propose Atlas Replay, an atlas-
based segmentation approach that uses prototypes to generate high-quality
segmentation masks through image registration that maintain consistency even
as the training distribution changes. We explore how our proposed method
performs compared to state-of-the-art CL methods in terms of knowledge
transferability across seven publicly available prostate segmentation
datasets. Prostate segmentation plays a vital role in diagnosing prostate
cancer, however, it poses challenges due to substantial anatomical variations,
benign structural differences in older age groups, and fluctuating acquisition
parameters. Our results show that Atlas Replay is both robust and generalizes
well to yet-unseen domains while being able to maintain knowledge, unlike end-
to-end segmentation methods. Our code base is available under
https://github.com/MECLabTUDA/Atlas-Replay.

##  1 Introduction

Continual learning (CL) plays a crucial role in safety-critical applications
of Deep Learning (DL), particularly in healthcare. In such domains, models
must continually adapt to data drift over time while maintaining high
performance on older data, even in cases where direct access to part of the
data is restricted for privacy reasons. The objective of CL is to train a
model that demonstrates high performance on sequentially arriving datasets,
despite the constrained timeframe during which the datasets are accessible.
Achieving this objective is challenging, as approaches tend to fall into one
of the extremes: either suffering from _catastrophic forgetting_ [19] by being
too plastic, or _unable to learn_ new tasks by being too rigid. Additionally,
some methods exhibit linear growth in training time and resource requirements
as the number of training tasks increases. CL models are further susceptible
to domain shifts over time [15], increasing the amount of catastrophic
forgetting even further. While various strategies have been proposed for CL,
they often fail to perform well with medical data [9], resulting in
segmentations that do not meet basic semantic standards. Therefore, striking a
delicate balance between preserving previous knowledge and maintaining the
necessary plasticity to learn new tasks is a key requirement when training in
a continual fashion [8, 16].

As we have domain knowledge on the geometry of the organ to be segmented, we
go back to the roots of medical image segmentation and contextualize atlas-
based segmentation [26], as an alternative to end-to-end CL segmentation
pipelines. In the context of CL, the access to structural information is key
for achieving accurate and semantically coherent predictions. In atlas-based
segmentation, registration is used to adapt the best-fitting labeled mask from
an atlas (i.e. a pre-defined set of reference images) to a patient scan. The
natural reliance of atlas-based approaches on geometrical aspects makes it
ideally applicable for CL setups as the structural information can be
leveraged to extract domain knowledge.

(a) Privacy concerns

(b) Rigidity/plasticity trade-off – No semantic coherence

(c) Dependent on capacity

Figure 1: Demonstration of different problems CL methods face.

CL for end-to-end models can be generally classified into one of the following
three categories: (1) (pseudo-) rehearsal, (2) regularization/distillation,
and (3) expansion. Figure 1 demonstrates the practical downsides of each
strategy.

The most successful CL methods follow some type of replay or rehearsal, which
involves storing samples to interleave them during later training [33, 37].
Long-term storage of patient scans, however, violates data protection
regulations [12]. Distillation methods such as _PLOP_ [11] do not store data
directly but rather distill knowledge from previously trained models. Yet
recent work [14] shows that the additional computational burden associated
with the pseudo-data generation hinders their use with high-dimensional
medical images. Regularization methods [6, 19] on the other hand work by
penalizing severe shifts from the previously-learned parameter space. These
approaches have lower resource requirements, but disappointing performance
across all tasks and merely allow a trade-off between rigidity and plasticity.
Models with high plasticity increase the ability to learn new information,
whereas rigid models maintain more knowledge from previous tasks and prevent
catastrophic forgetting. Additionally, models trained end-to-end in a
sequential manner with regularization or distillation approaches tend to
generate predictions with no semantic coherence. Expansion techniques maintain
stable performance across all tasks but grow the model size with the number of
tasks [15, 17]. Other works [29, 31, 32] propose architectural modifications
that are useful in certain settings but imply a high computational overhead.
Instead of relying on end-to-end CL, we propose a modular atlas-based method
for continual segmentation. _Atlas Replay_ leverages structural information to
maintain knowledge over time with the benefits of rehearsal while preserving
patient privacy and avoiding model growing.

Atlas Replay generates _prototypes_ built from a set of patient scans and the
_VoxelMorph_ registration framework [4] to perform registration between a
specific patient scan and a prototype. A prototype is a combination of
multiple patient scans – and corresponding segmentation masks – that disallows
the direct identification of a subject. Prototypes are registered to new
images to generate segmentation masks for new patients during deployment.

We introduce an approach to build an atlas of prototypes, irrespective of the
anatomy; and propose a CL method to perform atlas-based segmentation which
maintains model plasticity while preserving previous knowledge and outperforms
state-of-the-art (SOTA) end-to-end continual segmentation approaches. We can
consider using a stored atlas to generate segmentation masks as a form of
pseudo-rehearsal that maintains relevant information from previous examples
without storing actual patient images.

The contributions of this work are three-fold. Our proposed approach:

  * •

succesfully builds privacy-preserved prototypes and therefore being more
protective towards patient privacy compared to traditional replay-based
methods,

  * •

leverages prototypes and structural information to maintain knowledge over
time,

  * •

benefits from a rehearsal based approach and VoxelMorph to achieve stable
performances for continual image segmentation.

To validate our method, we investigate the problem of prostate segmentation in
T2-weighted Magnetic Resonance Images (MRI), which is an important step in the
diagnosis and treatment of prostate cancer [35]. Prostates have relatively
static shapes for which domain knowledge can be leveraged over time.
Variations in imaging protocol, such as the diminishing use of endorectal coil
over time [20], lead to domain shift. Such shifts clearly state the importance
of CL from a clinical perspective.

With the introduction of Atlas Replay, we pave the path for integrating atlas-
based methods into the realm of CL. We demonstrate that established
conventional approaches like atlas-based segmentation, which have fallen into
a certain neglect due to the current DL era, can be effectively utilized in
dynamic clinical setups. This success can be attributed to the use of
structural information, which showcases highly favorable outcomes in CL
scenarios.

##  2 Methodology

Traditional end-to-end segmentation methods in CL typically struggle to strike
a suitable balance between rigidity and plasticity, which presents a
significant drawback in this field. Such models often face a trade-off where
they either fail to acquire new knowledge by retaining information from
previously seen data, or overly prioritize recent cases, leading to
catastrophic forgetting and predictions lacking semantic coherence.
Recognizing the importance of structural information for achieving a proper
trade-off between rigidity and plasticity in CL setups, we turn to atlas-based
segmentation as a solution. Most atlas-based segmentation methods are based on
traditional (non-DL) techniques [1, 26]. We combine the advantages of DL with
the ability to leverage domain knowledge of atlas-based segmentation for CL.
An atlas $\Lambda$ is a manually labeled set of patient scans [34]. Atlas-
based segmentation is the modular process of registering an image from the
atlas to a new scan to directly generate an accurate segmentation by
transforming the respective mask in the same manner [26], thus leveraging
structural information. This approach has been successfully applied to
prostate cancer [3], heart regions [13], brain tissue [1] and aortic tissue
[18] for MRI scans.

#### Fundamentals

We start by introducing some key terminology: $\Omega\subset\mathbb{R}^{3}$
defines a 3D spatial domain. $\mathcal{T}_{p}\subset\Omega_{\mathcal{T}}$ is
referred to as dataset $p$ and consists of $m\\{(f^{i}_{m},f^{s}_{m})\\}$
pairs, where $f^{i}$ is a patient scan and $f^{s}$ the corresponding
segmentation mask. Stage $x$ in a continual setup defines the process of
training the model on dataset $x$ after it has been trained on all previous
$\\{1,\dots,x-1\\}$ tasks. $\Omega_{\mathcal{T}}$ is a set of datasets and
$\Omega_{\mathcal{P}}$, a set of prototypes. A prototype $\mathcal{P}_{k}$ is
a tuple $\mathcal{P}_{k}=(\mathcal{P}^{i}_{k},\mathcal{P}^{s}_{k})$, where
$\mathcal{P}^{i}_{k}$ is the scan and $\mathcal{P}^{s}_{k}$ the corresponding
segmentation mask. We define our set of prototypes $\Omega_{\mathcal{P}}$ as
privacy-preserving representations.

_VoxelMorph_ [4] is a popular framework for medical image registration. The
underlying architecture is a simple U-Net, to which additional convolutional
layers are added to generate a deformation field $\phi$. The network is
trained by penalizing the difference between the warped moving image using
$\phi$ and the target image. By using $\phi$ to alter the prototype
segmentation mask, one can utilize the framework for atlas-based segmentation
given a specific atlas $\Lambda$.

#### Prototype building

We create four distinct prototypes based on the coil type used during
acquisition. For each prototype $\mathcal{P}_{k}$, we extract $r$ random
samples from the associated dataset $\mathcal{T}_{p}$, a subset of
$\Omega_{\mathcal{T}}$. In this study, we set $r=7$, corresponding to the size
of the smallest training set. To evaluate the performance between prototypes,
we allocate three datasets from our data corpus to validate the inter-
prototype performance in Section 4.4.

First, $r$ random images are selected from the training dataset
$\mathcal{T}_{p}$. The first image $I$ represents the initial prototype
$\mathcal{P}_{k}$. For each following image $I^{*}$, we refine
$\mathcal{P}_{k}$ through rigid alignment using _SimpleITK_ [27, 38]. We then
compute the average of $\mathcal{P}_{k}$ and $I^{*}$ to update the prototype,
including the corresponding segmentation masks. The final $\mathcal{P}k$ is
the prototype for dataset $\mathcal{T}{p}$. Taking the average over multiple
scans hinders the recovery of patient-sensitive information during storing or
interleaving, as illustrated in Figure 2.

Figure 2: Illustration of prototype $\mathcal{P}_{A}$.

Figure 2 shows the difference between the final prototype $\mathcal{P}_{D}$
and three out of seven subject scans that were used to build the prototype.
The intensity distribution for every dataset is shown in Figure 3. We asses
the effectiveness of our prototypes in maintaining privacy by performing a
user study among senior radiologists, with more than 10 years of experience,
in Section 4.3.

Figure 3: Intensity histograms for every dataset.

#### Prototype registration

A patient scan $f^{i}$ represents the fixed image that is registered to the
best-fitting prototype $\mathcal{P}_{k}$ – the moving image – using
VoxelMorph. The deformation field $\phi$ from the network is then used to warp
the prototypes‘ segmentation mask: $\mathcal{P}_{k}^{s}\circ\phi$. As
illustrated in Figure 4, our method only stores prototypes that preserve key
information and therefore preserves patient privacy better than storing actual
subjects.

Figure 4: Comparison of sequential training on $n$ tasks against traditional
rehearsal and _Atlas Replay_ in terms of privacy preservation;
$\Omega_{\mathcal{P}}$ represents the stored prototypes.

#### Training continually using Atlas Replay

Algorithm 1 demonstrates a continual training setup in simplified pseudo-code
for our proposed Atlas Replay approach using the VoxelMorph (VxM) Framework.

Input: Datasets to train on
$\\{\mathcal{T}_{p}\\}_{p\leq\lvert\Omega_{\mathcal{T}}\rvert}$

Output: Trained model weights $\theta$

// Initialize $\mathcal{M}_{\theta}$

1 $\theta\leftarrow initializeModel()$

// Build $k=4$ prototypes

2 $\\{\mathcal{P}_{k}\\}\leftarrow
buildPrototypes\left(\\{\mathcal{T}_{p}\\}_{p\leq\lvert\Omega_{\mathcal{T}}\rvert},k=4\right)$

// Select prototype for $\mathcal{T}_{1}$

3 $\mathcal{P}\in\\{\mathcal{P}_{k}\\}$

// Train with VxM

4 $\theta\leftarrow\text{VxM}(\theta,\mathcal{T}_{1},\mathcal{P})$

5for  _$i\leftarrow 2$ to $\lvert\Omega_{\mathcal{T}}\rvert$_ do

6

// Select prototype for $\mathcal{T}_{i}$

7 $\mathcal{P}\in\\{\mathcal{P}_{k}\\}$

// Train with VxM

8 $\theta\leftarrow\text{VxM}(\theta,\mathcal{T}_{i},\mathcal{P})$

9 end for

Algorithm 1 Training using Atlas Replay

##  3 Experimental Setup

In this section, we briefly describe our corpus of publicly available datasets
and report essential aspects of our experimental setup.

#### Data corpus

Our prostate data corpus $\Omega_{\mathcal{T}}$ consists of seven publicly
available datasets [22, 5, 21, 23, 24, 2]. Table 1 shows the number of cases
in each dataset and the random 80:20 train/validation split. The splits along
with our code base is accessible under https://github.com/MECLabTUDA/Atlas-
Replay.

Dataset Task ID # Cases (train, val) Vendor Source RUNMC $\mathcal{T}_{1}$ 30
– (24, 6) Siemens [5] BMC $\mathcal{T}_{2}$ 30 – (24, 6) Philips HCRUDB
$\mathcal{T}_{3}$ 19 – (15, 4) Siemens [21] UCL $\mathcal{T}_{4}$ 13 – (10, 3)
Siemens [22] BIDMC $\mathcal{T}_{5}$ 12 – (9, 3) GE HK $\mathcal{T}_{6}$ 12 –
(9, 3) Siemens DecathProst $\mathcal{T}_{7}$ 32 – (25, 7) Unknown [2]

Table 1: Characteristics of our prostate data corpus; including the vendor of
the acquisition device.

#### Prototypes

To give the reader a proper understanding of our set of prototypes
$\Omega_{\mathcal{P}}$, these are illustrated in Figure 5.

(a) $\mathcal{P}_{A}$

(b) $\mathcal{P}_{B}$

(c) $\mathcal{P}_{C}$

(d) $\mathcal{P}_{D}$

Figure 5: Prototypes and their overlapping segmentation mask. For every
prototype, a random slice was selected.

#### Training setup

The VoxelMorph framework [4] with default optimizer and learning rate is used
for all registration experiments. To ensure a fair comparison between the
network’s performance, the same underlying U-Net architecture is used for end-
to-end segmentation. Registration models are trained for 250 epochs using a
composition of the Normalized Cross Correlation ($\mathcal{L}_{NCC}$), Cross-
Entropy (CE, $\mathcal{L}_{CE}$), and a gradient smoothing loss for the
deformation field ($\mathcal{L}_{smooth}$). The total registration loss is
calculated in the following way:

$\displaystyle\mathcal{L}_{reg}=\mathcal{L}_{NCC}+2\cdot\mathcal{L}_{CE}+\mathcal{L}_{smooth}$
(1)

Based on our ablation results from Section 4.5, we use the double-weighted CE
loss setup for all our VoxelMorph experiments. Segmentation networks are
trained for 250 epochs with $\mathcal{L}_{CE}$ only. All experiments were
carried out on a GeForce RTX 3090 GPU (24 GB).

#### Metrics

We report the mean Dice and standard deviation across the test images from all
tasks as well as average backward (BWT) and forward (FWT) transferability [25,
10]. BWT indicates the amount of maintained knowledge on test samples
$(f^{i}_{m},f^{s}_{m})\subset\mathcal{T}_{j}$ during training on different
stages $\\{\mathcal{T}_{p}\\}_{p\leq\lvert\Omega_{\mathcal{T}}\rvert}\;;\;j<p$
over time. FWT on the other hand measures the impact of the current training
stage $\\{\mathcal{T}_{p}\\}_{p\leq\lvert\Omega_{\mathcal{T}}\rvert}$ on test
data $(f^{i}_{m},f^{s}_{m})\subset\mathcal{T}_{j}\;;\;j>p$ from an untrained
stage.

Let $\mathcal{T}_{p}$ be a specific task:

$\displaystyle\text{BWT}\left(\mathcal{T}_{p}\right)$
$\displaystyle=\text{Dice}\left(\mathcal{M}_{\left[\mathcal{T}_{1},\dots,\mathcal{T}_{p},\dots,\mathcal{T}_{n}\right]},\mathcal{T}_{p}\right)$
$\displaystyle\>-\text{Dice}\left(\mathcal{M}_{\left[\mathcal{T}_{1},\dots,\mathcal{T}_{p}\right]},\mathcal{T}_{p}\right),$
(2)

where $\mathcal{M}_{\left[\mathcal{T}_{1},\dots,\mathcal{T}_{p}\right]}$ is a
network trained on stages
$\\{1,\dots,p\\}\leq\lvert\Omega_{\mathcal{T}}\rvert$ and
$\text{Dice}(\mathcal{M}_{\left[\mathcal{T}_{1},\dots,\mathcal{T}_{j}\right]},\mathcal{T}_{p})$
indicates the Sørensen–Dice coefficient from a network trained on stages
$\\{1,\dots,j\\}$ evaluated on dataset $p$. FWT is defined as:

$\displaystyle\text{FWT}\left(\mathcal{T}_{p}\right)$
$\displaystyle=\text{Dice}\left(\mathcal{M}_{\left[\mathcal{T}_{1},\dots,\mathcal{T}_{p-1}\right]},\mathcal{T}_{p}\right)$
$\displaystyle\>-\text{Dice}\left(\mathcal{M}_{\left[\mathcal{T}_{p}\right]},\mathcal{T}_{p}\right).$
(3)

FWT for the last model state as well as BWT for the first model state are not
defined. Models with high plasticity are able to learn new knowledge and
achieve higher FWT, whereas models that maintain most knowledge from previous
tasks obtain a higher BWT.

To validate our user study, we report sensitivity, specificity, precision
(positive predicted value), and the Matthews Correlation Coefficient (MCC)
[36, 7, 28]. We calculate the MCC to provide an overall assessment of the
classification performance with respect to both, true negative and true
positive rates.

#### Baselines

We compare Atlas Replay to end-to-end sequential training, the upper bound of
simple replay training, which requires the storage of actual patient scans,
and four popular CL methods: Elastic Weight Consolidation (EWC) [19],
Riemannian Walk (RWalk) [6], Incremental Learning Techniques (ILT) [30] using
distillation on the output (KD), intermediate (MSE), or both (KD, MSE) layers
and Bias Correction (BiC) [37]. A hyperparameter search is conducted for EWC
and RWalk, and the best settings are used of $\lambda=2.2$ (EWC),
$\alpha=0.9$, and $\lambda=1.7$ (RWalk). For ILT, we used the default
distillation parameter $\lambda_{D}=1$. Since BiC is a rehearsal-based method,
from each task, seven samples were interleaved.

##  4 Results

We present a comprehensive evaluation of various aspects of our work. We start
by analyzing the performance of Atlas Replay compared to sequential training,
EWC, RWalk, BiC, and rehearsal training. We then proceed with a qualitative
temporal evaluation in Section 4.2. Section 4.3 explores the effectiveness of
privacy-preserving prototypes through a user study conducted with senior
radiologists and computer scientists. We assess the inter-prototype
performance during inference in 4.4. Furthermore, we conduct a loss ablation
study in 4.5 and compare the U-Net’s end-to-end segmentation with our atlas-
based segmentation approach in 4.6. The results are based on the
hyperparameter search we conducted in Section 4.7. These evaluations provide
multiple insights into the performance, robustness, and effectiveness of the
proposed approach and prototype generation strategy.

###  4.1 Continual learning performance

In this section, we compare Atlas Replay to training a U-Net model
sequentially alongside the continual learning methods EWC, RWalk, ILT and BiC.
Further, we compare against the upper bound of rehearsal training (storing
seven samples from each task).

(a) Sequential

(b) EWC

(c) RWalk

(d) BiC

(e) Rehearsal

(f) Atlas Replay

Figure 6: CL performance as Dice; the larger the area the better.

Figure 6 shows the mean Dice scores of the final networks evaluated across all
seven datasets. EWC (6(b)) achieves more or less the same Dice performance
across the tasks as the simple sequential setup (6(a)) showing slight shifts
in BWT and FWT as seen in Table 2. RWalk, on the other hand, achieves a better
BWT, meaning it maintains more knowledge. However, this comes at the cost of a
lower mean Dice and FWT (6(c)). This indicates that even after tuning the
hyperparameters, end-to-end state-of-the-art CL methods under-perform for
medical image segmentation. BiC on the other hand preserves previous knowledge
the best (BWT), but at the cost of plasticity loss (FWT), Table 2. Atlas
Replay (6(f)) is the only method that maintains good performance over all seen
tasks. BWT is higher than for the rehearsal upper bound, as indicated in 6(d)
to 6(f), without compromising the plasticity loss (FWT) as observed by BiC,
Table 2. The different ILT versions on the other hand achieve slightly better
performances than EWC.

Method BWT $\uparrow$ [%] FWT $\uparrow$ [%] Sequential $-26.60\pm 17.88$
$-25.10\pm 17.44$ EWC $-27.48\pm 15.00$ $-29.65\pm 18.32$ RWalk $-22.79\pm
18.94$ $-51.48\pm 13.23$ ILT${}_{\text{KD}}$ $-27.86\pm 11.85$ $-26.79\pm
21.32$ ILT${}_{\text{MSE}}$ $-25.56\pm 14.90$ $-34.26\pm 18.94$
ILT${}_{\text{KD, MSE}}$ $-16.59\pm 12.77$ $-34.58\pm 17.73$ BiC
$\mathbf{-0.56\pm 7.42}$ $-29.01\pm 7.42$ Rehearsal $-6.60\pm 7.02$ $-24.23\pm
15.58$ Atlas Replay $-8.13\pm 7.68$ $\mathbf{-18.26\pm 15.22}$

Table 2: CL performance indicated by BWT and FWT along with standard
deviation.

With the default U-Net specifications of VoxelMorph, the average mean
segmentation Dice is about $61\%$. However, a significant improvement in our
proposed method can be clearly observed.

###  4.2 Qualitative temporal evaluation

To analyze the robustness of our proposed method qualitatively, we visualize
segmentation masks in Figure 7.

Figure 7: Temporal analysis for sequential, EWC, and Atlas Replay using Case
3, Slice 19 from $\mathcal{T}_{3}$.

Atlas Replay consistently produces coherent segmentation masks, irrespective
of the training stage. Sequential training and EWC on the other hand produce
low-quality segmentations until the network is trained on the particular stage
3 $\left(\\{\mathcal{T}_{p}\\}_{p\leq 3}\right)$. In particular, the low
performance on later tasks shows the effect of catastrophic forgetting, where
the network adapts too strongly to the later training data. Atlas Replay is
neither too rigid nor plastic, as is outputs robust predictions for data from
both early and later training stages.

The image shows the impact of selecting the correct prototype during inference
for the effectiveness of the registration method. Using a prototype that was
built from a non-coil dataset ($\mathcal{P}_{A}$) to perform registration with
a coil-acquired sample ($\mathcal{T}_{2}$) is expected to have lower
performance as the geometric shape of the prostate changes significantly
depending on the type of coil. Since the coil type is a central aspect of the
examination that can be easily recognized, it is to be expected that the user
applying the algorithm knows the type of coil that was used during acquisition
and can select the correct prototype accordingly. If this were not known, a
simple solution would be to register the desired sample across all prototypes
and then use the prototype leading to the best registration performance in
terms of NCC or MSE.

###  4.3 Effectiveness of privacy preserving prototypes

In order to evaluate the effectiveness of our privacy-preserving prototypes,
we conducted a user study involving two senior radiologists with more than 10
years of experience (RAD) and two computer scientists (CS). The study aimed to
assess the participants’ ability to correctly select the used subject from a
set of three patient scans, with only one scan being the correct match for the
shown prototype. This selection process was performed five times for each
prototype.

To measure the performance of the participants, we report four key metrics:
sensitivity, specificity, precision (positive predicted value), and Matthews
Correlation Coefficient (MCC). Sensitivity refers to the ability to correctly
identify the true positive cases, i.e., correctly selecting the matching
patient scan. Specificity measures the ability to correctly identify the true
negative cases, meaning to correctly exclude the non-matching patient scans.
Precision represents the proportion of correctly selected matching patient
scans out of the total selected matching scans. MCC on the other hand provides
an overall assessment of the classification performance, taking into account
both true positive and true negative rates. Table 3 shows the average results
of our user study for the radiologists and computer scientists including the
standard deviation and random chance.

Group Sensitivity $\uparrow$ [%] Specificity $\uparrow$ [%] Precision
$\uparrow$ [%] MCC $\uparrow$ [%] random – – $33.33\left(\frac{1}{3}\right)$ –
chance RAD $\mathbf{47.50\pm 7.50}$ $\mathbf{74.00\pm 4.00}$ $\mathbf{47.50\pm
7.50}$ $\mathbf{21.50\pm 14.50}$ CS $45.00\pm 10.00$ $73.00\pm 5.00$ $45.00\pm
10.00$ $18.00\pm 15.00$

Table 3: Results of our user study evaluation on privacy-preserving prototypes
compared to random chance. The evaluation included participation from both
radiologists and computer scientists.

Table 3 clearly shows that correctly selecting the used subject from a set of
three samples, along with the prototype, proved to be very challenging for
senior radiologists as well as computer scientists. The results illustrate
that possessing technical knowledge about the process of prototype building
has minimal impact on selecting the correct samples as the classification
rates are very similar to the ones from the radiologists. Such a difficulty in
correctly identifying the matching patient scan demonstrates the effectiveness
of our privacy-preserving prototype building approach. The average precision
achieved by the radiologists was $47.50\%$, indicating the complexity of the
task and the privacy preservation capabilities of our prototypes.

###  4.4 Inter-prototype performance

Figure 8: Combination of scatter and box plot showing performance distribution
and change in Dice for every task based on the used prototype during inference
with the final model trained sequentially on all datasets.

To assess the inter-prototype performance, we trained sequentially on all
seven tasks using Atlas Replay and utilized the final model to analyze the
influence of prototypes during inference (Figure 8). Specifically, for each
sample from the test sets we performed registration using all four prototypes.

Figure 8 illustrates the distribution of Dice for every task based on the
utilized prototype during inference with the sequentially trained model. This
visualization effectively shows how the choice of prototype during the
registration process influences the performance observed in inference across
different tasks.

###  4.5 Loss ablation study

To determine the optimal setup for atlas-based segmentation, we conduct an
ablation study where we modify the traditional registration loss proposed in
the VoxelMorph paper [4]. This study allows us to identify the most effective
VoxelMorph configuration. The networks are trained on the joint dataset
$\bigcup_{p=1}^{\lvert\Omega_{\mathcal{T}}\rvert}\mathcal{T}_{p}$ and the
evaluation is performed across all seven datasets
$\\{\mathcal{T}_{p}\\}_{p\leq\lvert\Omega_{\mathcal{T}}\rvert}$.

Ablation $\left(\mathcal{L}_{reg},\ref{align:loss}\right)$ Dice
$\uparrow{}\pm{}$ $\sigma\downarrow$ [%] $\mathcal{T}_{1}$ $\mathcal{T}_{2}$
$\mathcal{T}_{3}$ $\mathcal{T}_{4}$ $\mathcal{T}_{5}$ $\mathcal{T}_{6}$
$\mathcal{T}_{7}$
$\mathcal{L}_{NCC}+0\cdot\mathcal{L}_{CE}+\mathcal{L}_{smooth}$ $37.03\pm
17.05$ $62.63\pm 1.97$ $66.56\pm 6.06$ $69.78\pm 2.08$ $58.48\pm 6.94$
$63.73\pm 8.74$ $38.84\pm 11.86$
$\mathcal{L}_{NCC}+1\cdot\mathcal{L}_{CE}+\mathcal{L}_{smooth}$ $56.01\pm
15.05$ $68.95\pm 3.46$ $70.83\pm 3.36$ $73.80\pm 1.27$ $\mathbf{60.18\pm
8.94}$ $63.28\pm 11.71$ $57.89\pm 10.83$
$\mathcal{L}_{NCC}+2\cdot\mathcal{L}_{CE}+\mathcal{L}_{smooth}$
$\mathbf{70.55\pm 6.32}$ $\mathbf{73.29\pm 3.65}$ $\mathbf{71.39\pm 1.94}$
$\mathbf{79.84\pm 3.49}$ $60.07\pm 4.88$ $\mathbf{65.36\pm 8.89}$
$\mathbf{69.47\pm 6.70}$ $\mathcal{L}_{CE}$ (U-Net joint) $79.62\pm 2.50$
$77.01\pm 11.10$ $78.95\pm 6.61$ $79.72\pm 4.13$ $83.49\pm 1.97$ $82.45\pm
1.98$ $80.34\pm 1.94$

Table 4: VoxelMorph ablations trained on joint prostate data with different
loss variations (top) and U-Net joint results (bottom).

Table 4 shows that weighting the Cross-Entropy loss twice leads to the best
segmentation performance across all tasks. A network trained with no
segmentation loss term, i.e. only for registration, achieves a limited
contribution in terms of segmentation performance. Given this insight, we use
the double weighted Cross-Entropy loss setup for all our VoxelMorph related
experiments as shown in Equation 1.

###  4.6 End-to-end segmentation vs. atlas-based segmentation

To better assess the generalizability of the models, we train U-Net and Atlas
Replay networks for every dataset $\mathcal{T}_{p}\subset\Omega_{\mathcal{T}}$
and validated them across all datasets.

Figure 9: Scatter plot showing the performance distribution in Dice when
trained on one dataset and evaluated across all seven datasets from the
prostate corpus $\Omega_{\mathcal{T}}$; trend line is based on mean Dice
performance across all evaluated samples.

Figure 9 illustrates the difference between an end-to-end segmentation
approach (U-Net) and Atlas Replay in terms of generalizability and
segmentation performance. Using the proposed registration-based approach, we
obtain robust results regardless of what dataset is used for training. The
U-Net performance on the other hand is clearly dependent on the training
dataset, as indicated by the strong distribution shifts in the figure, i.e.
lack of generalizability. $\mathcal{T}_{5}$, for instance, shows a significant
performance deterioration. Besides an increase in generalizability, the
performance of the registration-based method is also higher, as shown by the
mean Dice trend line, which is consistently above that of the U-Net. These
results show the increased versatility and robustness of Atlas Replay even
before continual learning is performed.

###  4.7 Hyperparameter search

For every used end-to-end CL method, a hyperparameter search was performed
using four different settings. The setup with the highest mean performance in
terms of Dice, BWT, and FWT together was selected for our experiments. Table 5
shows the results for each hyperparameter setting. The highest values for a
method and the best parameter settings are marked in bold.

Method Prostate Fixed params Tuned param Dice $\uparrow$ [%] BWT $\uparrow$
[%] FWT $\uparrow$ [%] EWC – $\lambda=0.4$ $42.44\pm 20.67$ $-28.47\pm 12.35$
$-28.87\pm 19.48$ $\lambda=1.1$ $44.83\pm 20.86$ $-29.28\pm 14.82$
$\mathbf{-27.47\pm 16.29}$ $\lambda=1.7$ $44.03\pm 20.49$ $-30.96\pm 15.33$
$-30.59\pm 16.66$ $\mathbf{\lambda=2.2}$ $\mathbf{46.06\pm 18.46}$
$\mathbf{-27.48\pm 15.00}$ $-29.65\pm 18.32$ RWalk $\alpha=0.9,$ update $=20$
$\lambda=0.4$ $29.94\pm 23.59$ $-29.01\pm 12.84$ $\mathbf{-46.69\pm 13.29}$
$\lambda=1.1$ $27.80\pm 20.07$ $-25.91\pm 17.21$ $-50.06\pm 6.757$
$\mathbf{\lambda=1.7}$ $\mathbf{30.58\pm 21.38}$ $-22.79\pm 18.94$ $-51.48\pm
13.23$ $\lambda=2.2$ $30.48\pm 23.26$ $\mathbf{-21.43\pm 10.50}$ $-54.31\pm
11.07$

Table 5: Results of the hyperparameter search considered for the two used end-
to-end CL methods; mean Dice, BWT and FWT over all tasks including standard
deviation [in %]; highest values and best parameter setup are marked in bold.

The results clearly show the addressed trade-off problem from Figure 1(b).
Depending on how the hyperparameter(s) are set, the trade-off between
plasticity and knowledge preservation varies. The network either performs well
in terms of maintaining knowledge over time – increased $\lambda$ for EWC – or
the results are very similar no matter how the hyperparameter is set –
$\lambda$ for RWalk.

##  5 Conclusion

We introduce Atlas Replay, a robust atlas-based segmentation technique for
continuous training in clinical environments with data drift. We evaluate our
approach on seven prostate segmentation scenarios and show that it outperforms
state-of-the-art continual learning methods by maintaining knowledge from
early stages without compromising model plasticity. _Atlas Replay_ is more
generalizable than U-Net models even when trained statically with only data
from one site, with a performance difference of $15\%$ on average. We
additionally introduced a prototype-building method with initial privacy
preservation from a human perspective that can be leveraged to maintain
structural information over time resulting in a proper balance between
rigidity and plasticity for CL setups. Future work should address more
stringent privacy preservation for prototypes and a self-growing atlas
technique in combination with proper prototype adjustments over time to
further push the potential of registration for segmentation in terms of CL. By
releasing our code base along with instructions and pre-trained networks, we
hope to amplify and inspire CL research beyond end-to-end models that
incorporates atlas-based segmentation for medical settings.

##  6 Acknowledgements

This work was (partially) supported by RACOON (NUM), under BMBF grant number
01KX2021 and EVA-KI, under BMG grant number ZMVI1- 2520DAT03A.

Supplementary Material

##  1 End-to-end vs. atlas-based segmentation

Table 6 provides the Dice scores with standard deviation for every trained
baseline evaluated across all tasks which were used to create Figure 9 from
the main manuscript.

Baselines Dice $\uparrow{}\pm{}$ $\sigma\downarrow$ [%] $\mathcal{T}_{1}$
$\mathcal{T}_{2}$ $\mathcal{T}_{3}$ $\mathcal{T}_{4}$ $\mathcal{T}_{5}$
$\mathcal{T}_{6}$ $\mathcal{T}_{7}$ U-Net $\mathcal{T}_{1}$ $\mathbf{75.38\pm
5.96}$ $37.82\pm 15.47$ $17.04\pm 4.83$ $42.02\pm 23.75$ $6.66\pm 6.86$
$50.86\pm 1.78$ $76.14\pm 5.66$ $\mathcal{T}_{2}$ $36.36\pm 19.55$
$\mathbf{70.25\pm 6.30}$ $25.81\pm 9.26$ $52.21\pm 15.53$ $50.30\pm 8.01$
$36.78\pm 7.98$ $37.01\pm 10.70$ $\mathcal{T}_{3}$ $21.12\pm 17.64$ $53.55\pm
6.23$ $\mathbf{78.50\pm 5.04}$ $55.44\pm 2.65$ $37.54\pm 6.98$ $50.12\pm 5.39$
$18.94\pm 9.90$ $\mathcal{T}_{4}$ $46.48\pm 17.13$ $65.69\pm 6.29$ $60.23\pm
2.08$ $\mathbf{66.44\pm 11.89}$ $49.19\pm 13.37$ $65.73\pm 10.23$ $44.92\pm
11.52$ $\mathcal{T}_{5}$ $9.61\pm 7.19$ $48.44\pm 8.02$ $29.82\pm 5.24$
$40.26\pm 18.72$ $\mathbf{59.55\pm 4.97}$ $38.21\pm 4.85$ $7.87\pm 7.35$
$\mathcal{T}_{6}$ $22.05\pm 19.25$ $54.02\pm 10.05$ $58.41\pm 4.21$ $59.37\pm
4.41$ $52.10\pm 1.99$ $\mathbf{77.42\pm 3.27}$ $22.13\pm 17.11$
$\mathcal{T}_{7}$ $73.52\pm 7.57$ $37.38\pm 9.17$ $20.28\pm 3.27$ $32.83\pm
22.86$ $10.75\pm 7.59$ $45.11\pm 3.84$ $\mathbf{76.12\pm 3.83}$ Atlas Replay
$\mathcal{T}_{1}$ $\mathbf{70.50\pm 4.63}$ $56.02\pm 4.78$ $72.87\pm 6.30$
$73.79\pm 1.55$ $34.81\pm 22.70$ $63.69\pm 9.80$ $73.02\pm 3.63$
$\mathcal{T}_{2}$ $44.27\pm 18.47$ $\mathbf{72.30\pm 6.05}$ $73.92\pm 3.50$
$75.16\pm 0.91$ $53.73\pm 13.33$ $66.62\pm 7.50$ $43.99\pm 9.74$
$\mathcal{T}_{3}$ $34.30\pm 7.69$ $34.69\pm 6.41$ $\mathbf{73.11\pm 4.42}$
$71.33\pm 2.28$ $23.81\pm 11.07$ $44.28\pm 8.05$ $31.97\pm 7.77$
$\mathcal{T}_{4}$ $23.20\pm 9.93$ $23.61\pm 7.13$ $63.36\pm 7.56$
$\mathbf{78.40\pm 0.72}$ $20.57\pm 7.71$ $37.11\pm 4.82$ $19.55\pm 5.34$
$\mathcal{T}_{5}$ $32.79\pm 15.13$ $53.30\pm 2.42$ $72.87\pm 5.67$ $74.55\pm
2.59$ $\mathbf{60.31\pm 3.02}$ $58.50\pm 9.19$ $30.31\pm 8.77$
$\mathcal{T}_{6}$ $34.66\pm 15.19$ $44.05\pm 6.22$ $68.91\pm 3.17$ $71.48\pm
4.08$ $39.16\pm 15.10$ $\mathbf{64.97\pm 13.40}$ $33.47\pm 8.32$
$\mathcal{T}_{7}$ $69.02\pm 5.74$ $53.63\pm 8.64$ $74.94\pm 5.43$ $73.69\pm
1.35$ $36.15\pm 16.09$ $59.41\pm 10.48$ $\mathbf{67.73\pm 5.20}$

Table 6: Results for all baseline networks trained on every task individually
and evaluated across all tasks; Bold values indicate the performance of the
baseline on the validation set of the task it has been trained on.

## References

  * [1] Paul Aljabar, Rolf A Heckemann, Alexander Hammers, Joseph V Hajnal, and Daniel Rueckert.  Multi-atlas based segmentation of brain images: atlas selection and its effect on accuracy.  Neuroimage, 46(3):726–738, 2009.
  * [2] Michela Antonelli, Annika Reinke, Spyridon Bakas, Keyvan Farahani, Bennett A Landman, Geert Litjens, Bjoern Menze, Olaf Ronneberger, Ronald M Summers, Bram van Ginneken, et al.  The medical segmentation decathlon.  arXiv preprint arXiv:2106.05735, 2021.
  * [3] Takahiro Aoyama, Hidetoshi Shimizu, Tomoki Kitagawa, Kazushi Yokoi, Yutaro Koide, Hiroyuki Tachibana, Kojiro Suzuki, and Takeshi Kodaira.  Comparison of atlas-based auto-segmentation accuracy for radiotherapy in prostate cancer.  Physics and Imaging in Radiation Oncology, 19:126–130, 2021.
  * [4] Guha Balakrishnan, Amy Zhao, Mert R Sabuncu, John Guttag, and Adrian V Dalca.  Voxelmorph: a learning framework for deformable medical image registration.  IEEE transactions on medical imaging, 38(8):1788–1800, 2019.
  * [5] Nicholas Bloch, Anant Madabhushi, Henkjan Huisman, et al.  Nci-isbi 2013 challenge: Automated segmentation of prostate structures.  Medical image analysis, 2015.
  * [6] Arslan Chaudhry, Puneet K Dokania, Thalaiyasingam Ajanthan, and Philip HS Torr.  Riemannian walk for incremental learning: Understanding forgetting and intransigence.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 532–547, 2018.
  * [7] Davide Chicco and Giuseppe Jurman.  The advantages of the matthews correlation coefficient (mcc) over f1 score and accuracy in binary classification evaluation.  BMC genomics, 21:1–13, 2020.
  * [8] Matthias De Lange, Rahaf Aljundi, Marc Masana, Sarah Parisot, Xu Jia, Aleš Leonardis, Gregory Slabaugh, and Tinne Tuytelaars.  A continual learning survey: Defying forgetting in classification tasks.  IEEE transactions on pattern analysis and machine intelligence, 44(7):3366–3385, 2021.
  * [9] Mohammad Mahdi Derakhshani, Ivona Najdenkoska, Tom van Sonsbeek, Xiantong Zhen, Dwarikanath Mahapatra, Marcel Worring, and Cees GM Snoek.  Lifelonger: A benchmark for continual disease classification.  -, 2022.
  * [10] Natalia Díaz-Rodríguez, Vincenzo Lomonaco, David Filliat, and Davide Maltoni.  Don’t forget, there is more than forgetting: new metrics for continual learning.  In Workshop on Continual Learning, NeurIPS 2018 (Neural Information Processing Systems, 2018.
  * [11] Arthur Douillard, Yifu Chen, Arnaud Dapogny, and Matthieu Cord.  Plop: Learning without forgetting for continual semantic segmentation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4040–4050, 2021.
  * [12] European Commission.  Regulation (EU) 2016/679 of the European Parliament and of the Council of 27 April 2016 on the protection of natural persons with regard to the processing of personal data and on the free movement of such data, and repealing Directive 95/46/EC (General Data Protection Regulation) (Text with EEA relevance), 2016.
  * [13] Tarun Kanti Ghosh, Md Kamrul Hasan, Shidhartho Roy, Md Ashraful Alam, Eklas Hossain, and Mohiuddin Ahmad.  Multi-class probabilistic atlas-based whole heart segmentation method in cardiac ct and mri.  IEEE Access, 9:66948–66964, 2021.
  * [14] Camila Gonzalez, Amin Ranem, Daniel Pinto dos Santos, Ahmed Othman, and Anirban Mukhopadhyay.  Lifelong nnunet: a framework for standardized medical continual learning.  -, 2022.
  * [15] Camila González, Amin Ranem, Ahmed Othman, and Anirban Mukhopadhyay.  Task-agnostic continual hippocampus segmentation for smooth population shifts.  In MICCAI Workshop on Domain Adaptation and Representation Transfer, pages 108–118. Springer, 2022.
  * [16] Raia Hadsell, Dushyant Rao, Andrei A Rusu, and Razvan Pascanu.  Embracing change: Continual learning in deep neural networks.  Trends in cognitive sciences, 24(12):1028–1040, 2020.
  * [17] Ching-Yi Hung, Cheng-Hao Tu, Cheng-En Wu, Chien-Hung Chen, Yi-Ming Chan, and Chu-Song Chen.  Compacting, picking and growing for unforgetting continual learning.  Advances in Neural Information Processing Systems, 32, 2019.
  * [18] Ivana Isgum, Marius Staring, Annemarieke Rutten, Mathias Prokop, Max A Viergever, and Bram Van Ginneken.  Multi-atlas-based segmentation with local decision fusion—application to cardiac and aortic segmentation in ct scans.  IEEE transactions on medical imaging, 28(7):1000–1010, 2009.
  * [19] James Kirkpatrick, Razvan Pascanu, Neil Rabinowitz, Joel Veness, Guillaume Desjardins, Andrei A Rusu, Kieran Milan, John Quan, Tiago Ramalho, Agnieszka Grabska-Barwinska, et al.  Overcoming catastrophic forgetting in neural networks.  Proceedings of the national academy of sciences, 114(13):3521–3526, 2017.
  * [20] Grace Lee, Aytekin Oto, and Mihai Giurcanu.  Prostate mri: Is endorectal coil necessary?—a review.  Life, 12(4):569, 2022.
  * [21] Guillaume Lemaître, Robert Martí, Jordi Freixenet, Joan C Vilanova, Paul M Walker, and Fabrice Meriaudeau.  Computer-aided detection and diagnosis for prostate cancer based on mono and multi-parametric mri: a review.  Computers in biology and medicine, 60:8–31, 2015.
  * [22] Geert Litjens, Robert Toth, Wendy Van De Ven, Caroline Hoeks, Sjoerd Kerkstra, Bram van Ginneken, Graham Vincent, Gwenael Guillard, Neil Birbeck, Jindang Zhang, et al.  Evaluation of prostate segmentation algorithms for mri: the promise12 challenge.  Medical image analysis, 18(2):359–373, 2014.
  * [23] Quande Liu, Qi Dou, and Pheng Ann Heng.  Shape-aware meta-learning for generalizing prostate mri segmentation to unseen domains.  In International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI), 2020.
  * [24] Quande Liu, Qi Dou, Lequan Yu, and Pheng Ann Heng.  Ms-net: Multi-site network for improving prostate segmentation with heterogeneous mri data.  IEEE Transactions on Medical Imaging, 2020.
  * [25] David Lopez-Paz and Marc’Aurelio Ranzato.  Gradient episodic memory for continual learning.  Advances in neural information processing systems, 30, 2017.
  * [26] Maria Lorenzo-Valdés, Gerardo I Sanchez-Ortiz, Raad Mohiaddin, and Daniel Rueckert.  Atlas-based segmentation and tracking of 3d cardiac mr images using non-rigid registration.  In International conference on medical image computing and computer-assisted intervention, pages 642–650. Springer, 2002.
  * [27] Bradley C Lowekamp, David T Chen, Luis Ibáñez, and Daniel Blezek.  The design of simpleitk.  Frontiers in neuroinformatics, 7:45, 2013.
  * [28] Brian W Matthews.  Comparison of the predicted and observed secondary structure of t4 phage lysozyme.  Biochimica et Biophysica Acta (BBA)-Protein Structure, 405(2):442–451, 1975.
  * [29] Marius Memmel, Camila Gonzalez, and Anirban Mukhopadhyay.  Adversarial continual learning for multi-domain hippocampal segmentation.  In Domain Adaptation and Representation Transfer, and Affordable Healthcare and AI for Resource Diverse Global Health, pages 35–45. Springer, 2021.
  * [30] Umberto Michieli and Pietro Zanuttigh.  Incremental learning techniques for semantic segmentation.  In Proceedings of the IEEE/CVF international conference on computer vision workshops, pages 0–0, 2019.
  * [31] Sinan Özgün, Anne-Marie Rickmann, Abhijit Guha Roy, and Christian Wachinger.  Importance driven continual learning for segmentation across domains.  In International Workshop on Machine Learning in Medical Imaging, pages 423–433. Springer, 2020.
  * [32] Amin Ranem, Camila González, and Anirban Mukhopadhyay.  Continual hippocampus segmentation with transformers.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 3711–3720, 2022.
  * [33] Sylvestre-Alvise Rebuffi, Alexander Kolesnikov, Georg Sperl, and Christoph H Lampert.  icarl: Incremental classifier and representation learning.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, pages 2001–2010, 2017.
  * [34] Torsten Rohlfing, Robert Brandt, Randolf Menzel, Daniel B Russakoff, and Calvin R Maurer.  Quo vadis, atlas-based segmentation?  In Handbook of biomedical image analysis, pages 435–486. Springer, 2005.
  * [35] Clare MC Tempany, Peter R Carroll, Michael S Leapman, W Robert Lee, and Diane MF Savarese.  The role of magnetic resonance imaging in prostate cancer.  UpToDate. Waltham (MA): UpToDate, 2018.
  * [36] Robert Trevethan.  Sensitivity, specificity, and predictive values: foundations, pliabilities, and pitfalls in research and practice.  Frontiers in public health, 5:307, 2017.
  * [37] Yue Wu, Yinpeng Chen, Lijuan Wang, Yuancheng Ye, Zicheng Liu, Yandong Guo, and Yun Fu.  Large scale incremental learning.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 374–382, 2019.
  * [38] Ziv Yaniv, Bradley C Lowekamp, Hans J Johnson, and Richard Beare.  Simpleitk image-analysis notebooks: a collaborative environment for education and reproducible research.  Journal of digital imaging, 31(3):290–303, 2018.


# A model for cooperative scientific research inspired by the ant colony
algorithm

Zhuoran He  School of Physics, Huazhong University of Science and Technology,
Wuhan 430074, China  Tingtao Zhou  Division of Engineering and Applied
Sciences, California Institute of Technology, Pasadena, CA 91125, USA

###### Abstract

Modern scientific research has become largely a cooperative activity in the
Internet age. We build a simulation model to understand the population-level
creativity based on the heuristic ant colony algorithm. Each researcher has
two heuristic parameters characterizing the goodness of his own judgments and
his trust on literature. In a population with all kinds of researchers, we
find that as the problem scale increases, the contributor distribution
significantly shifts from the independent regime of relying on one’s own
judgments to the cooperative regime of more closely following the literature.
The distribution also changes with the stage of the research problem and the
computing power available. Our work provides some preliminary understanding
and guidance for the dynamical process of cooperative scientific research in
various disciplines.

##  I Introduction

Cooperative scientific research is a new trend in the science community
nowadays due to the growth of number of researchers [1, 2, 3], the faster
propagation of knowledge through the Internet [4, 5, 6, 7] and the many new
interdisciplinary research topics [8, 9], etc. Research groups ranging from a
few scientists to international institutions can study related problems and
build upon each other’s works. In the early pioneering days, the activity of
scientific research was the solitary work of a few geniuses of the world and
the spirits of independent thinking and skepticism were highly valued. In
modern days, we are seeing more and more scientific achievements made by the
progressive efforts of many researchers. This paradigm shift accompanies the
development of complexity science itself [8, 9, 10]. Scientists in the
Internet age work like a highly cooperative ant colony connected by pheromone,
i.e., research publications, and exhibit population-level creativity which
requires modeling to understand and optimize.

Previous studies on scientific research and collective intelligence have
discussed various aspects of this topic including the citation system [11, 12,
13, 14], evaluation and funding system [15, 16], game theory competition and
cooperation [17, 18, 19, 20], complex networks [21, 22], team size and
composition management [23, 24, 25], and so on [26, 27, 28, 29, 30, 31]. In
this paper, we build a simplified model inspired by the ant colony
optimization (ACO) algorithm [32, 33, 34] to study the dynamical process of
cooperative scientific research by computer simulations. Our ant colony model
can obtain the optimal research styles for various types of scientific
problems, e.g., simple (elemental) v.s. complex, new v.s. old (long-standing),
etc., and study the influence of computing power and different survival rules
on selecting researchers for the community.

We suppose that in the ant community, each scientific problem they study is a
randomly generated traveling salesman problem (TSP) [35] with $N$ vertices,
where $N$ controls the complexity of the problem. A researcher’s effort on
such a problem is modeled as making small decisions step by step to connect
the vertices and find a plausible path. He will then pass on the knowledge by
publications, i.e., leaving pheromone on the edges visited. The shorter the
total path, the more pheromone will be assigned. Every decision made is
governed by two heuristic parameters: $\alpha$ characterizing the researcher’s
trust on published literature, i.e., the pheromone left on all edges, and
$\beta$ characterizing his trust on the greedy local distance measure, i.e.,
the researcher’s own sense of direction. The procedure is iterated as
generations of researchers attempt for better solutions. Finally, the
accumulated pheromone concentrates on the shortest TSP path found by the
community, which represents the currently best answer known to the scientific
problem.

Two essential ingredients of our ant colony model are the NP-hardness of TSP
and the pheromone mechanism in ACO. Since TSP is NP-hard, it is easy to
evaluate and compare path lengths and exclude the longer path as ‘wrong’, but
difficult to know if the shorter path is indeed shortest, which is similar to
open questions in science that satisfy the falsifiability criterion. The
pheromone is a population-level information sharing mechanism that enables
researchers to work out difficult scientific problems cooperatively. Our ant
colony model develops the ACO in that we have improved the pheromone update
rules and allow the heuristic parameters $\alpha,\beta$ to differ individually
and change by evolution. We can then study the equilibrium distributions of
$\alpha,\beta$ given different problem scales $N$ and different numbers of ACO
iterations that distinguish between new and old problems. The influence of
computing power will be modeled by introducing the Hamiltonian cycle speedup
[36] that mimics the role of computers.

##  II Ant colony model

###  II.1 The core ACO algorithm

In the ACO algorithm, each ant with two heuristic parameters $\alpha,\beta$
tries to find a TSP path individually. The ant has its own memory of the set
of vertices $S$ that has been visited and has access to the community-shared
information $d_{ij}$, the distance between vertices $i,j$ and $\tau_{ij}$, the
amount of pheromone on the undirectioned edge $i-j$. The ant picks a random
vertex to start the trip. Then each step from vertex $i$ to vertex $j$ is
determined by the transition probability

$P_{i\rightarrow
j}\propto\,\begin{cases}\;(0.01+\tau_{ij})^{\alpha}\,/\,d_{ij}^{\beta},&j\not\in
S,\\\ \;0,&j\in S.\end{cases}$ (1)

The probabilities $P_{i\rightarrow j}$ are normalized for all $j\not\in S$ to
determine the next stop $j$. Then vertex $j$ is added to set $S$ so it will
not be repeatedly visited. Eq. (1) describes the basic rules of the
heuristically biased self-avoiding walk (SAW) [37, 38] of ants in the original
ACO algorithm. We have added a small background value of $0.01$ to $\tau_{ij}$
in Eq. (1) so that the ants do not get oversensitive to small amounts of
pheromone.

We have also made improvements in the pheromone update rules. After
$N_{\mathrm{ants}}=50$ ants have finished their TSP paths, we pick the best
$p(t)$ percent and allow these winning ants to leave pheromone over their TSP
paths bidirectionally. The amount of pheromone on each edge is inversely
proportional to the total path length and proportional to a linearly
decreasing weight of the path ranking. Long-distance steps on the TSP path
gets extra penalties. The pheromone on all edges then evaporates by $p(t)$
percent and the above procedure iterates while the percentage $p(t)$ gradually
decreases from 50% to 8% (4 ants) over the iterations.

These improvements mean that initially the ant colony is very eager to
accumulate pheromone and later, the update rules get tighter as the best-known
path of the ant colony becomes nearly optimal. But any time, an ant who beats
the best-known path always immediately becomes the leader of the top 4 ants
and leaves the most pheromone to the whole ant colony. With the improved
pheromone update rules, the pheromone becomes a more useful guide to the ants
and better resembles the literature publication system in academia. More
details of the model can be found in our Matlab code provided upon reasonable
request.

###  II.2 Evolution of heuristic parameters

In the original ACO algorithm, the parameters $\alpha,\beta$ were set manually
as hyper-parameters and applied to all ants. In our model, we allow
$\alpha,\beta$ to take different values for different ants and evolve the
distribution $P_{\alpha,\beta}$ by training the ant colony with randomly
generated TSP graphs. During the solution of one graph, $P_{\alpha,\beta}$ is
kept unchanged and the heuristic parameters of the ants who found shorter TSP
paths than the best-known path are recorded. These ants are called
contributors and their $\alpha,\beta$ values are used for evolving
$P_{\alpha,\beta}$ according to

$P_{\alpha,\beta}^{\mathrm{(new)}}=\frac{n_{c}}{M}P_{\alpha,\beta}^{(c)}+\left(1-\frac{n_{c}}{M}\right)P_{\alpha,\beta}^{\mathrm{(old)}}.$
(2)

Here $n_{c}$ is the number of contributors recorded during one graph and
$M=4000$ is the total pool of ants out of which the $N_{\mathrm{ants}}=50$
ants are sampled in each ACO iteration. When equilibrium is reached, every ant
in the colony is equally likely to become a contributor. More favorable
$(\alpha,\beta)$ values will attract more ants and less favorable values will
be adopted by fewer ants.

We then consider a more sophisticated situation where the trained distribution
$P_{\alpha,\beta}(t)$ can depend on problem stage $t$. To do this, we record
for each contributor not only its $\alpha,\beta$ values, but also the number
of ACO iterations $t$ performed when its contribution is made. We can then
compare at equilibrium the distributions $P_{\alpha,\beta}(t)$ suitable for
different problem stages $t$.

###  II.3 Hamiltonian cycle speedup

The Hamiltonian cycle speedup [36] is often used in conjunction with ACO to
speed up its convergence. In the core ACO algorithm, at every vertex $i$, the
ant heuristically picks a next vertex $j$ based on $P_{i\rightarrow j}$, which
mimics human intuition. The TSP path obtained is called a Hamiltonian cycle in
graph theory, e.g.,

$1\rightarrow(i-1)\rightarrow i\rightarrow\cdots\rightarrow
j\rightarrow(j+1)\rightarrow N\rightarrow 1.$ (3)

The Hamiltonian cycle speedup plays the role of a computer exhaustively
checking human errors. It enumerates all segments
$i\rightarrow\cdots\rightarrow j$ of the Hamiltonian cycle in Eq. (3) and
checks if the cycle length can be made shorter by reversing the segment into
$j\rightarrow\cdots\rightarrow i$, which is true if and only if
$d_{i-1,i}+d_{j,j+1}>d_{i-1,j}+d_{i,j+1}$. The exhaustive check continues
until no such improvements are possible, which is a necessary but not
sufficient condition for the optimality of the TSP path. We examine the
influence of the Hamiltonian cycle speedup on the distribution
$P_{\alpha,\beta}$ in Sec. III C and use the simpler model described in Secs.
II A–II B elsewhere.

##  III Results and discussions

###  III.1 Effect of problem scale

Figure 1: Histograms of heuristic parameters $\alpha$ (blue) and $\beta$ (red)
and contour plots of their joint distributions (insets) for TSP graphs with
$N=10$, 20, 50, 100, 200, and $500$ vertices. The “+” sign marks the mode
peak.

We first examine how the equilibrium distribution $P_{\alpha,\beta}$ changes
with problem scale $N$, i.e., the number of vertices in the TSP graph. The
vertices are randomly sampled from a uniform distribution in a 2D unit square
region. We have tried other distributions (Gaussian, triangular) and other
region shapes (rectangle, circle) and have found qualitatively the same
results.

As shown in Figs. 1a–1f, the $\alpha$ peak significantly shifts to larger
values as the problem scale $N$ increases. This indicates that when faced with
more difficult problems, the research community has to rely more on previous
works for guidance to find sophisticated better solutions, and random trials
ignoring literature is not as efficient. The $\beta$ parameter governs the
researcher’s goodness of local distance measure or sense of direction. For
small problems, the optimal distribution relies heavily on high $\beta$
values. For larger problems, the weights of high $\beta$ reach a plateau and
the joint distribution $P_{\alpha,\beta}$ develops a positive correlation
between $\alpha$ and $\beta$, suggesting that the successful research style is
a combination of the $\alpha$ and $\beta$ heuristics. Such researchers always
keeps up-to-date knowledge of the currently best solution known by the
community ($\alpha$ heuristics) and quickly identifies where potential
improvements are possible ($\beta$ heuristics) around the community-found
path.

###  III.2 Time-dependent distribution

Figure 2: (a) The $P_{\alpha,\beta}(t)$ distribution with $t=1$–$5$, $6$–$30$,
$31$–$100$ and $101$–$1000$ scatter plotted in blue to red. (b) Error curves
of $P_{\alpha,\beta}(t)$ compared with $P_{\alpha,\beta}$ and other
variations. Panels (c)–(d) compare the $P_{\alpha,\beta}$ distributions with
and without the Hamiltonian cycle speedup. Number of vertices $N=100$ in all
$4$ subplots.

Some difficult problems can persist for years or decades as researchers come
and leave. In an ideal situation, researchers switching from problems to
problems specialize to both the appropriate scale of complexity and the stage
of problem conducive to their own research styles ($\alpha,\beta$ values). We
therefore consider $P_{\alpha,\beta}(t)$ with $t$ being the number of ACO
iterations for fixed problem scale $N=100$. We train $P_{\alpha,\beta}(t)$ to
equilibrium and plot the results in Fig. 2a in 4 colors corresponding to 4
problem stages: newly proposed ($t=1$–$5$, blue), early ($t=6$–$30$, green),
intermediate ($t=31$–$100$, yellow), and late ($t=101$–$1000$, red) periods.

When a problem is newly proposed, the contributors (blue) generally have high
$\beta$ values. Since there are not many publications to read yet, researchers
with low $\beta$ will move randomly between the vertices and obtain TSP paths
of order $\mathcal{O}(N)$, while those with high $\beta$ will always greedily
choose the closest vertex to move to. The greedy solution can be estimated to
be

$\mathcal{O}\left(\frac{1}{\sqrt{N}}+\frac{1}{\sqrt{N-1}}+\cdots+1\right)=\mathcal{O}(\sqrt{N}),$
(4)

which is much better than a random self-avoiding walk $\mathcal{O}(N)$.
Therefore, all contributors of a newly proposed problem tend to have high
$\beta$ values. After the greedy solution has been found, the early-stage
contributors (green) constitute the upper part of the time-independent
distribution $P_{\alpha,\beta}$ in Fig. 1d. The intermediate (yellow) and
late-stage (red) contributors then scan down to the lower part of
$P_{\alpha,\beta}$ and finally concentrate into a red blob below the mode peak
of $P_{\alpha,\beta}$.

The error curves of the time-dependent $P_{\alpha,\beta}(t)$ and time-
independent $P_{\alpha,\beta}$ distributions are compared in Fig. 2b. We use
the relative error $\varepsilon(t)=L(t)/L_{\mathrm{opt}}-1$ averaged over 500
graphs to evaluate the goodness of a given research condition, where
$L_{\mathrm{opt}}$ is the optimal path length obtained from the open-source
exact TSP solver Concorde [39, 40]. Initially, the greedy solution of
$P_{\alpha,\beta}(t)$ (yellow dashed line in Fig. 2b) has an advantage over
$P_{\alpha,\beta}$ (blue line in Fig. 2b), which does not last for very long.
The blue and yellow curves nearly coincide when the problem reaches
intermediate to late periods. The residue error at $t=10^{3}$ remains $\sim
1.2\%$, which is likely to result from the path dependence effect [41], i.e.,
the ant colony gets trapped to a local minimum found by previous works. If we
have two or more independent research communities (green & purple lines in
Fig. 2b), which is realized by running the ant colony code multiple times and
keeping the smallest TSP length of the trials at every iteration step $t$, the
relative error $\varepsilon(t)$ has a statistically significant reduction.

###  III.3 Effects of computing power

We then move on to discuss the effects of more computing power, which is
mimicked by introducing the Hamiltonian cycle speedup described in Sec. II C.
When individual researchers have computers that help them do exhaustive trials
and verifications, our results indicate that the selectivity effects on both
the literature parameter $\alpha$ and the intuition parameter $\beta$ of the
contributors are significantly reduced. As is shown in Figs. 2c-2d, both the
$\alpha$ peak and the $\beta$ plateaus are made lower by introducing the
Hamiltonian cycle speedup. This means that computers are a chance equalizer
which diversifies the heuristic parameter distributions of the contributors.
Also, the $\alpha$ peak slightly shifts to smaller values, which is due to the
reduction of effective problem hardness when computers become available. In
terms of relative error (red line in Fig. 2b), the introduction of Hamiltonian
cycle speedup significantly speeds up the convergence to the optimal TSP path.
The residue error at $t=10^{3}$ iterations $\sim 0.4\%$ is made much smaller
than the blue curve but still nonzero, which suggests that the path dependence
effect of ant-colony research cannot be completely eliminated even with more
computing power available to each researcher individually.

###  III.4 Other non-ideal situations

Figure 3: (a) The equilibrium contributor distribution under the survival
rules of $t\leq 50$ (blue) and $t>50$ (red). Dashed contours are those of the
standard $N=100$ distribution in Fig. 1d. (b) The error curves
$\varepsilon(t)$ of different survival rules. Inset shows normally trained
contributor distribution in terms of improvement percentage v.s. problem
stage.

We often see in academia that researchers are faced with tight and pressing
survival rules, most of which are achievement-based. We find in our model that
sometimes these rules can be counter-productive to the science community. An
important reason why this happens is that such rules would encourage
researchers to focus on new or early-stage problems, leaving the late-stage
problems simply “outdated” rather than solved.

We simulate such a situation and results are shown in Fig. 3. Suppose a
problem is interesting to the ant colony for only $t\leq 50$ iterations, after
which the problem becomes old and out of attention. By training the ant colony
using TSP graphs with $N=100$ vertices under such hasty rules, the equilibrium
distribution $P_{\alpha\beta}$ is given by the blue dots in Fig. 3a.
Conversely, if every graph is solved up to $t=10^{3}$ iterations but only
those contributors after $t>50$ are recorded to update $P_{\alpha\beta}$, the
ant colony will be trained into the distribution of the red dots. The green
contours in Fig. 3a show the normally trained distribution where all
contributors are recorded to update $P_{\alpha\beta}$. Since achievement-based
survival rules pick out those contributors with big improvements of TSP
lengths, which, according to the inset of Fig. 3b, tend to be early-stage
contributors, the distribution $P_{\alpha\beta}$ shifts to the blue side as a
result.

We then plot in Fig. 3b the error curves $\varepsilon(t)$ of different
$P_{\alpha\beta}$ distributions averaged over 500 graphs. The blue
distribution has short-term benefits but long-term costs. A “hasty” ant colony
adapted to early-stage problems would lack those ants with heuristic
parameters conducive to making breakthroughs on long-standing problems and
therefore become inefficient as problems approach late stages. In reality, the
combined effect of making the researcher community both inefficient and not
interested in solving long-standing problems could be even worse, which can be
mimicked by reducing $N_{\mathrm{ants}}=50\,e^{-(t-50)/200}$ after $t>50$
(blue dashed line in Fig. 3b). The residue error at $t=10^{3}$ reaches $\sim
2\%$. More interestingly, the normally trained $N=100$ distribution (green
line in Fig. 3b) can be outperformed by the red distribution (red line in Fig.
3b) in the long run, which suggests the importance of giving more weights to
the late contributors.

##  IV Conclusion

We have established an ant-colony research model which enables us to
understand the dynamical process of cooperative scientific research in various
disciplines. Based on our model, we have made several interesting findings.
First, as the problem scale increases, the contributors tend to have more
cooperative heuristic parameters than those of simpler problems. Therefore,
the cooperative mode of scientific research is a consequence of complexity
science itself. Second, different problem stages will require different
research styles. In the beginning, simple intuitive thinking can help lay down
the general framework. Later, improvements become harder and require deeper
thinking and more trials and errors. Third, the introduction of computers or
any other advanced technology that enables exhaustive trials and verifications
can give the human researcher more freedom, diversify the contributor
population and make the scientific results more accurate and objective.

In addition to demonstrating the power of cooperative scientific research, our
model can also simulate non-ideal situations and identify how things might go
wrong. First is path dependence. As scientists build upon each other’s works,
there is inevitably some degree of path dependence. Parallel development of
several independent communities, technological methods, or schools of thoughts
can be better than having one unified community stuck with pre-established
ideas and paradigms. Second is hasty research. Putting pressure on
productivity or individual achievements can lead to hasty research. It is
important to give more credits to late contributors and solvers of long-
standing problems for the long-term progress of science.

## Acknowledgments

We thank the helpful discussions with colleagues in Department of physics at
Caltech. T. Z. acknowledges support by the Cecil and Sally Drinkward
Fellowship.

## References

  * Bloom  _et al._ [2020] N. Bloom, C. I. Jones, J. Van Reenen, and M. Webb, Are ideas getting harder to find?, American Economic Review 110, 1104 (2020).
  * Larson  _et al._ [2014] R. C. Larson, N. Ghaffarzadegan, and Y. Xue, Too many phd graduates or too few academic job openings: The basic reproductive number r0 in academia, Systems research and behavioral science 31, 745 (2014).
  * Gibbons  _et al._ [2018] M. Gibbons, B. L. Yoder, and J. Roy, Engineering by the numbers, in  _American Society for Engineering Education_ (2009, 2012, 2018) p. 37.
  * Ginsparg [2011] P. Ginsparg, Arxiv at 20, Nature 476, 145 (2011).
  * Goodrum  _et al._ [2001] A. A. Goodrum, K. W. McCain, S. Lawrence, and C. L. Giles, Scholarly publishing in the internet age: a citation analysis of computer science literature, Information Processing & Management 37, 661 (2001).
  * Holmberg and Thelwall [2014] K. Holmberg and M. Thelwall, Disciplinary differences in twitter scholarly communication, Scientometrics 101, 1027 (2014).
  * Hurd [2000] J. M. Hurd, The transformation of scientific communication: A model for 2020, Journal of the American society for information science 51, 1279 (2000).
  * Yin  _et al._ [2017] G.-Z. Yin, W.-B. Zhang, and S. Z. Cheng, Giant molecules: where chemistry, physics, and bio-science meet, Science China Chemistry 60, 338 (2017).
  * Bazzani  _et al._ [2009] A. Bazzani, B. Giorgini, and S. Rambaldi, Traffic and crowd dynamics: The physics of the city, in  _Encyclopedia of Complexity and Systems Science_, edited by R. A. Meyers (Springer New York, New York, NY, 2009) pp. 9411–9429.
  * Anderson [1972] P. W. Anderson, More is different, Science 177, 393 (1972).
  * Garfield [1955] E. Garfield, Citation indexes for science, Science 122, 108 (1955).
  * Leydesdorff [1998] L. Leydesdorff, Theories of citation?, Scientometrics 43, 5 (1998).
  * Hirsch [2007] J. E. Hirsch, Does the h index have predictive power?, Proceedings of the National Academy of Sciences 104, 19193 (2007).
  * Fong and Wilhite [2017] E. A. Fong and A. W. Wilhite, Authorship and citation manipulation in academic research, PloS one 12, e0187394 (2017).
  * Hicks [2012] D. Hicks, Performance-based university research funding systems, Research policy 41, 251 (2012).
  * Muscio  _et al._ [2013] A. Muscio, D. Quaglione, and G. Vallanti, Does government funding complement or substitute private research funding to universities?, Research Policy 42, 63 (2013).
  * Lim and Wittek [2018] I. S. Lim and P. Wittek, Satisfied-defect, unsatisfied-cooperate: An evolutionary dynamics of cooperation led by aspiration, Physical Review E 98, 062113 (2018).
  * [18] L. Tiokhin, Competition for priority harms the reliability of science but reforms can help, Nature Human Behaviour 20.
  * Sonubi  _et al._ [2016] A. Sonubi, A. Arcagni, S. Stefani, and M. Ausloos, Effects of competition and cooperation interaction between agents on networks in the presence of a market capacity, Physical Review E 94, 022303 (2016).
  * Axelrod [1997] R. Axelrod,  _The complexity of cooperation: Agent-based models of competition and collaboration_ , Vol. 3 (Princeton university press, 1997).
  * Barrat  _et al._ [2004] A. Barrat, M. Barthelemy, R. Pastor-Satorras, and A. Vespignani, The architecture of complex weighted networks, Proceedings of the national academy of sciences 101, 3747 (2004).
  * Newman [2001] M. E. J. Newman, Clustering and preferential attachment in growing networks, Phys. Rev. E 64, 025102(R) (2001).
  * Milojević [2014] S. Milojević, Principles of scientific research team formation and evolution, Proceedings of the National Academy of Sciences 111, 3984 (2014).
  * Massey  _et al._ [2006] C. Massey, F. Alpass, R. Flett, K. Lewis, S. Morriss, and F. Sligo, Crossing fields: The case of a multi-disciplinary research team, Qualitative Research 6, 131 (2006).
  * Wu  _et al._ [2019] L. Wu, D. Wang, and J. A. Evans, Large teams develop and small teams disrupt science and technology, Nature 566, 378 (2019).
  * Goldman and Shaked [1991] A. I. Goldman and M. Shaked, An economic model of scientific activity and truth acquisition, Philosophical Studies 63, 31 (1991).
  * Kealey and Nelson [1996] T. Kealey and R. R. Nelson,  _The economic laws of scientific research_ (Macmillan London, 1996).
  * Wu  _et al._ [2015] Z.-X. Wu, Z. Rong, and H.-X. Yang, Impact of heterogeneous activity and community structure on the evolutionary success of cooperators in social networks, Physical Review E 91, 012802 (2015).
  * Wakeling and Bak [2001] J. Wakeling and P. Bak, Intelligent systems in the context of surrounding environment, Physical Review E 64, 051920 (2001).
  * Karamched  _et al._ [2020] B. Karamched, M. Stickler, W. Ott, B. Lindner, Z. P. Kilpatrick, and K. Josić, Heterogeneity improves speed and accuracy in social networks, Physical Review Letters 125, 218302 (2020).
  * Durve  _et al._ [2020] M. Durve, L. Piro, M. Cencini, L. Biferale, and A. Celani, Collective olfactory search in a turbulent environment, Physical Review E 102, 012402 (2020).
  * Dorigo [1992] M. Dorigo, Optimization, learning and natural algorithms, PhD Thesis, Politecnico di Milano (1992).
  * Dorigo  _et al._ [1991] M. Dorigo, A. Colorni, and V. Maniezzo, Distributed optimization by ant colonies (1991).
  * Dorigo  _et al._ [1996] M. Dorigo, V. Maniezzo, and A. Colorni, Ant system: optimization by a colony of cooperating agents, IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 26, 29 (1996).
  * Flood [1956] M. M. Flood, The traveling-salesman problem, Operations research 4, 61 (1956).
  * Croes [1958] G. A. Croes, A method for solving traveling-salesman problems, Operations Research 6, 791 (1958).
  * Slade [2019] G. Slade, Self-avoiding walk on the complete graph (2019), arXiv:1904.11149 [math.PR] .
  * Li [2020] Z. Li, Positive speed self-avoiding walks on graphs with more than one end, Journal of Combinatorial Theory, Series A 175, 105257 (2020).
  * Applegate  _et al._ [2006] D. Applegate, R. Bixby, V. Chvatal, and W. Cook, Concorde tsp solver (2006).
  * Applegate  _et al._ [2009] D. L. Applegate, R. E. Bixby, V. Chvátal, W. Cook, D. G. Espinoza, M. Goycoolea, and K. Helsgaun, Certification of an optimal tsp tour through 85,900 cities, Operations Research Letters 37, 11 (2009).
  * David [1985] P. David, Clio and the economics of qwerty, The American Economic Review 75, 332 (1985).