Further author information: (Send correspondence to Abhinav K. Jha)
Abhinav K. Jha<EMAIL_ADDRESS>This manuscript has been accepted to
SPIE Medical Imaging, February 19-23, 2023. Please use the following reference
when citing the manuscript.

Yu, Z., Rahman, M. A., Abbey, C. K., Siegel, B. A., and Jha, A. K.,
“Development and task-based evaluation of a scatter-window projection and deep
learning-based transmission-less attenuation compensation method for
myocardial perfusion SPECT”, Proc. SPIE Medical Imaging, 2023.

# Development and task-based evaluation of a scatter-window projection and
deep learning-based transmission-less attenuation compensation method for
myocardial perfusion SPECT

Zitong Yu  Department of Biomedical Engineering, Washington University in St.
Louis, St. Louis, USA  Md Ashequr Rahman  Department of Biomedical
Engineering, Washington University in St. Louis, St. Louis, USA  Craig K.
Abbey  Department of Psychological & Brain Sciences, University of California
Santa Barbara, Santa Barbara, USA  Barry A. Siegel  Mallinckrodt Institute of
Radiology, Washington University in St. Louis, St. Louis, USA  Abhinav K. Jha

###### Abstract

Attenuation compensation (AC) is beneficial for visual interpretation tasks in
single-photon emission computed tomography (SPECT) myocardial perfusion
imaging (MPI). However, traditional AC methods require the availability of a
transmission scan, most often a CT scan. This approach has the disadvantage of
increased radiation dose, increased scanner costs, and the possibility of
inaccurate diagnosis in cases of misregistration between the SPECT and CT
images. Further, many SPECT systems do not include a CT component. To address
these issues, we developed a Scatter-window projection and deep Learning-based
AC (SLAC) method to perform AC without a separate transmission scan. To
investigate the clinical efficacy of this method, we then objectively
evaluated the performance of this method on the clinical task of detecting
perfusion defects on MPI in a retrospective study with anonymized clinical
SPECT/CT stress MPI images. The proposed method was compared with CT-based AC
(CTAC) and no-AC (NAC) methods. Our results showed that the SLAC method
yielded an almost overlapping receiver operating characteristic (ROC) plot and
a similar area under the ROC (AUC) to the CTAC method on this task. These
results demonstrate the capability of the SLAC method for transmission-less AC
in SPECT and motivate further clinical evaluation.

###### keywords:

Attenuation compensation, Deep learning, SPECT reconstruction, Task-based
evaluation

##  1 INTRODUCTION

Attenuation of photons is a major image-degrading effect that adversely
impacts image quality in single-photon emission computed tomography (SPECT).
Multiple studies have shown that attenuation compensation (AC) is beneficial
for clinical interpretations of SPECT myocardial perfusion images.[1, 2]
Conventional AC methods typically require an attenuation map, now most
commonly obtained from a separate CT scan. However, these CT-based AC (CTAC)
methods have multiple disadvantages, such as increased radiation dose, higher
scanner costs, and possible misregistration between the SPECT and CT images
potentially leading to inaccurate diagnosis.[3, 4, 5, 6] Further, many SPECT
systems often do not have a CT component. For example, SPECT systems in
smaller community hospitals and physician offices, as well as mobile SPECT
systems facilitating use in remote locations are often SPECT-only. The
emerging solid-state-detector-based SPECT systems, which provide higher
sensitivity, energy, temporal, and spatial resolution compared to conventional
SPECT systems, often do not have CT imaging capability either.[7, 8] For these
reasons, there is an important need to develop transmission-less AC (Tx-less
AC) methods for SPECT.

Given this high significance, multiple Tx-less AC methods have been proposed,
including methods that use SPECT emission data to estimate attenuation maps[9,
10, 11] and methods that operate on the iterative inversion of the forward
mathematical models of SPECT systems.[12, 13] More recently, deep learning
(DL)-based methods have shown significant promise for Tx-less AC.[14, 15, 16,
17, 18, 19] Shi et al. recently reported promising performance of a
conditional generative adversarial network for Tx-less AC for myocardial
perfusion SPECT (MPS).[18] Chen et al. developed strategies for generating
attenuation maps using emission data for dedicated cardiac SPECT with small
field-of-view and found that their strategies outperformed the method that
directly predicts AC images from non-attenuation-corrected images.[19] While
the performance of these methods is promising, these DL-based methods have
typically been evaluated using figures of merit (FoM) that measure the
fidelity between the images reconstructed using the DL-based approach with a
reference standard, which is typically the image reconstructed with the CT-
based AC method. Medical images are acquired for specific clinical tasks.
Thus, clinical translation of these Tx-less AC methods requires that they be
evaluated in reference to the clinical task.[20, 21, 22, 23] However, studies
have shown that evaluation using fidelity-based FoMs may not correlate with
performance on clinical tasks in myocardial perfusion imaging (MPI).[24, 25,
26] Thus, it is crucial to evaluate these methods on the specific clinical
tasks for which the images are acquired.

We have shown that scatter-window data in SPECT contains information to
estimate the attenuation distribution.[27] Based on this premise, we had
proposed a DL-based Tx-less AC method for SPECT.[28] In this paper, we advance
upon this idea to propose a Scatter-window projection and DL-based Tx-less AC
(SLAC) method for myocardial perfusion SPECT (MPS) that uses only SPECT
emission data in photopeak and scatter windows. We objectively evaluate the
method on the clinical task of myocardial defect detection in a retrospective
study.

##  2 METHODS

###  2.1 Proposed method

The overall framework of the SLAC method is shown in Fig. 1. The probability
of scatter at a certain location is proportional to the attenuation
distribution at that location. It is expected that a reconstruction of the
scatter-window projection would show the contrast between regions with
different attenuation coefficients. A previous study has shown a promising
performance of a DL-based method that estimates the attenuation map from a
reconstruction of scatter-window projection.[18] Thus, the scatter-window
projection was reconstructed using an ordered-subsets expectation maximization
(OSEM)-based approach, yielding an initial estimate of attenuation map.[29,
30] Then, in this study, we used a DL-based technique to segment the initial
estimate of attenuation maps. U-Net-based approaches have shown promise in
biomedical image segmentation problems.[31, 32, 33] Thus, we used a U-Net-
based approach, namely, multi-channel input and multi-encoder U-Net (McEUN),
to segment the initial estimate of attenuation maps. The McEUN was trained to
segment the initial estimate of attenuation maps into six regions, including
skin and subcutaneous adipose, muscles and organs, lungs, bones, patient
table, and background. The McEUN mainly consists of two components: an encoder
with muti-channel input and an assembly of six decoders. To stabilize the
network training and leverage salient regions, skip connections with attention
gate (AG) were implemented between the output of layers in the encoder and
each decoder.[34] Dropout was applied to prevent overfitting.[35] The network
was designed to input the whole 3-D image and maximize the amount of
information learned in a global sense. The network was trained to minimize the
cross entropy between estimated and true segmentation.

Figure 1: The overall framework of the SLAC method.

###  2.2 Evaluation

We evaluated the SLAC method on the clinical task of detecting cardiac
perfusion defects in a retrospective Institutional Review Board-approved
evaluation study with anonymized stress SPECT MPI data. We compared the
performance of our method with activity maps reconstructed using a CT-based AC
(CTAC) method and to images obtained without AC, referred to as no-AC (NAC)
approach. We followed the recently proposed Recommendations for EvaLuation of
AI in NuClear-medicinE (RELAINCE) guidelines to lend rigor to our
evaluation[23]. The description of the evaluation study consists of four
parts, including data collection and curation, network training and method
implementation, process to extract task-specific information, and figures of
merit used.

####  2.2.1 Data collection and curation

The dataset used in this study consisted of N = 648 anonymized clinical
SPECT/CT stress MPI studies scanned between January 2016 and July 2021, with
SPECT projection data and CT images along with clinical reports. As per the
clinical reports, we categorized patients diagnosed with normal rest and
stress myocardial perfusion function as healthy patients, while patients
diagnosed with ischemia in a left ventricular wall as diseased patients. MPI
scans were acquired on a GE Discovery NM/CT 670 system after the injection of
${}^{\text{99m}}$Tc-tetrofosmin. SPECT emission data were collected in
photopeak (126-154 keV) and scatter windows (114-126 keV). CT images were
acquired at 120 kVp on a GE Optima CT 540 system integrated with the GE
Discovery NM/CT 670. To avoid misalignment between CT and SPECT scans, CT
images were registered to the SPECT space using MIM Maestro (MIM Software Inc,
Cleveland, OH). CT-defined attenuation maps were calculated from the CT scans
using a bi-linear model.[36]

For evaluating the SLAC method on the clinical task of detecting cardiac
defects, knowledge of the existence and location of the defects was needed.
The clinical records have limitations in providing this information. To
address this issue, we implemented a strategy that introduces synthetic
cardiac defects in healthy patient images.[37] We designed 27 types of
clinically realistic defects with three radial extents, three severities, and
at three locations. A summary of defect types is shown in Table 1.

Table 1: Defect parameters Parameter |
---|---
Extent | 30, 60, and 90 degrees around the left ventricular (LV) wall
Severity | 10%, 25%, and 50% less activity than the normal myocardium
Location | Anterior, inferior, and lateral LV walls

The whole dataset was divided into the training dataset (N = 508) and the
testing dataset (N = 140). Because we need the ground truth of defects in the
test dataset for our evaluation study, the testing dataset only consisted of
patients that were diagnosed as healthy according to the clinical records, and
synthetic defects were introduced in 71 of these 140 healthy test patients,
referred to as defect-present samples. The remaining 69 healthy test patients
were referred to as defect-absent samples. We generated 27×71 = 1917 defect-
present samples in both photopeak and scatter energy windows. We also
generated 27×69 = 1863 defect-absent samples, although they were identical if
from the same patient.

As mentioned in Sec. 2.1, the scatter-window projection was reconstructed
using an OSEM-based reconstruction method without AC, yielding the initial
estimate of the attenuation map. Also, the photopeak-window projection was
reconstructed using the same strategy, yielding an initial estimate of the
activity map. The CT- based attenuation maps for network training were
segmented into skin and subcutaneous adipose, muscles and organs, lungs,
bones, patient table, and background, using a Markov random field-based
method.[38] The average attenuation coefficients of each region were
calculated and served as the predefined attenuation coefficients.

####  2.2.2 Network training and method implementation

A total of N = 508 samples were used for the network training. The kernel
weights of the McEUN were initialized using the Glorot normal initializer.[39]
Biases were initialized to a constant of 0.03. The McEUN was trained to
minimize a weighted cross entropy loss between predicted and CT-based
segmentations using Adam optimizer.[40] We optimized the weight parameters to
yield the best segmentation performance. Five-fold cross-validation was
implemented to prevent overfitting. The training and validation were performed
using Keras 2.2.4 on two TITAN RTX GPUs with 23 GB memory each.

As mentioned in Sec. 2.2.1, there were a total of 27×140 = 3780 samples in the
test dataset. The trained McEUN yielded segmented masks of initial estimates
of attenuation maps. Predefined attenuation coefficients were assigned to each
region, yielding the final estimate of attenuation maps. Next, the photopeak-
window projections were reconstructed using an OSEM-based reconstruction
method, which accounts for attenuation and collimator-detector response.[30]
The final estimates of attenuation maps were used for AC. The reconstructed
activity images had a size of 64×64×64 with a voxel size of 0.68 cm. Following
the clinical protocols, the reconstructed activity images were reoriented into
short-axis slices and filtered by a Butterworth filter with an order of 5 and
cutoff frequency 0.44 cm-1. CTAC-based images and NAC-based images were
obtained using the same OSEM-based reconstruction approach and post-processing
procedures as used in the SLAC method but with different AC approaches.

####  2.2.3 Process to extract task-specific information

We objectively evaluated the performance of the SLAC method on the task of
detecting myocardial perfusion defects in an observer study. While ideally,
such evaluation should be performed with human observers, this is time-
consuming and tedious. Model observers provide an easy-to-use in silico
approach to perform such evaluation and identify methods for evaluation with
human observers. Thus, multiple studies use model observers to evaluate
imaging systems and methods.[21, 41, 42, 43] In our evaluation study, the
location of the defect was chosen to be in the inferior left ventricular wall
in the test dataset. Then, we had a patient population where the defect extent
and severity were varying, but the location of the defect was the same in the
entire population in the test dataset. Therefore, there were 9×71 = 639
defect-present samples and 9×69 = 621 defect-absent samples in the evaluation
study. Previous studies have shown that the channelized Hotelling observer
(CHO) with rotationally symmetric frequency channels can emulate human-
observer performance on the task of detecting perfusion defects from MPS
images in this setting.[44, 45] Thus, in this study, we used this model
observer. We extracted a 32×32 region from the middle 2-D slice of the short-
axis images that had the defect centroid at the center and applied the CHO to
this region to yield test statistics. These statistics were calculated for
each defect-present and -absent image in the test set using a leave-one-out
strategy.

####  2.2.4 Figures of merit

The test statistics generated were compared to a threshold to classify the
image into the defect-present or defect-absent class. By varying the
threshold, the ROC curve was plotted[46, 47] using LABROC4 program.[48] The
AUC measures the performance of methods on the task of defect detection. A
higher AUC indicates better performance. We calculated the AUC with 95%
confidence intervals (CIs) for SLAC, CTAC, and NAC methods.

To assess the performance of our method using visual fidelity-based criterion,
we calculated root mean-square-error (RMSE) and structural similarity index
(SSIM) with CIs between images obtained using the SLAC and the CTAC method, as
well as between images obtained using the NAC and the CTAC method.

##  3 RESULTS

###  3.1 Comparing to other AC methods

Fig. 2 shows ROC curves obtained by SLAC, CTAC, and NAC methods, along with
the corresponding 95% CIs. The ROC curve obtained by the SLAC method almost
overlapped those obtained by the CTAC method and outperformed the NAC method.
Fig. 3 shows the AUC values with CIs obtained by three methods. We observed
that the AUC obtained by SLAC was similar to that obtained by the CTAC method
and significantly outperformed ($p<0.05$) the NAC method.

Figure 2: ROC curves obtained by CTAC, SLAC, and NAC methods. Shadows indicate
95% confidence intervals.

Figure 3: AUC obtained by CTAC, SLAC, and NAC methods.

###  3.2 Representative examples

Fig. 4 shows examples of SPECT images and corresponding attenuation maps
estimated by SLAC, compared with those from CTAC. We found that the
attenuation maps obtained using SLAC method were close to those obtained from
CT images. Further, the short-axis SPECT images obtained using SLAC methods
were similar to those obtained using CTAC methods. Fig. 4.b shows the same
defect was shown in CTAC and SLAC-based images.

Figure 4: Examples of SPECT images and attenuation maps obtained by SLAC and
CTAC method: (a) a defect-absent example, (b) a defect-present example, where
the yellow arrow indicates the detect.

###  3.3 Evaluation using fidelity-based figures of merit

Table 2 shows RMSE and SSIM between images obtained using the SLAC and the
CTAC method, as well as between the NAC and the CTAC method. We found that the
SLAC method significantly outperformed the NAC method based on these metrics.

Table 2: Evaluation using fidelity-based FoMs Method | RMSE (95% CIs) | SSIM (95% CIs)
---|---|---
SLAC | 0.013 (0.011, 0.014) | 0.96 (0.95,0.97)
NAC | 0.036 (0.032,0.040) | 0.66 (0.65,0.68)

##  4 DISCUSSIONS AND CONCLUSIONS

This study represents a step along the pathway for developing a method for AC
in SPECT without requiring a transmission scan. Our results demonstrate that
the proposed SLAC method did not just yield activity images that were visually
similar to those yielded by the CTAC method (as evaluated using the RMSE and
SSIM metrics) (Table 2 and Fig. 4), but also yielded a similar performance to
the CTAC method on the clinical task of detecting perfusion defects from MPI
images (Fig. 3). More specifically, the proposed method yielded an almost
overlapping ROC and similar AUC to the CTAC method and significantly
outperformed the NAC method. Further, evaluation using fidelity-based FoMs
showed that the SLAC method significantly outperformed the NAC method (Table
2). These results demonstrate that the SLAC method has the capability of
performing AC using only emission data for SPECT.

Our evaluation results motivate evaluation of the SLAC method with data from
different scanners and from different institutions to assess for
generalizability. Another area of future research is advancing the method to
use list-mode data instead of sinogram data to reconstruct the attenuation
map. Previous studies have shown that processing data in list-mode format in
SPECT can yield improved performance on clinical tasks compared to processing
data in binned format.[49, 50, 51, 52, 53] More specifically, it has been
quantitatively shown that list-mode data contains more information compared to
binned data for the task of estimating attenuation coefficients.[27] Thus,
advancing the method to process data in list-mode format directly may lead to
even more improved performance.

There are some limitations to this study. First, the performance of the SLAC
method is dependent on the quality of training data, including the
segmentation of CT images. We used a Markov random field-based method to
segment the CT images. However, there are other methods, including DL-based
methods, which may yield segmented low-dose CT images with better segmentation
performance.[54, 55] Integrating a more advanced segmentation method is a
future direction for improving the performance of the SLAC method. Second, we
used a model observer study to evaluate the defect detection performance of
the SLAC method. While the observer we used has been shown to mimic human
observer performance, ideally this study should be conducted using human
observers, such as experienced radiologists. The results from the model
observer study motivate the evaluation of the method using human observers.
Furthermore, the test dataset only consisted of defects located in the
inferior LV wall. A future research direction is to incorporate defects at
multiple locations in the test dataset while performing the observer study.

In conclusion, a scatter window projection and deep learning-based
transmission-less AC method for SPECT yielded similar performance compared to
the standard CTAC method, as evaluated on the clinically relevant task of
detecting myocardial perfusion defects with a model observer in a
retrospective evaluation study.

###### Acknowledgements.

This work was partly supported by National Institute of Biomedical Imaging and
Bioengineering of National Institute of Health (NIH) under grant number
R21-EB024647, R01-EB031051, and R01EB031962. Support is also acknowledged from
the NVIDIA GPU grant. We also thank the Washington University Center for High
Performance Computing for providing computational resources. The center is
partially funded by NIH grants 1S10RR022984-01A1 and 1S10OD018091-01. We also
thank Dr. Craig K. Abbey for helpful discussions on the design of the observer
study.

## References

  * [1] Hutton, B. F., Buvat, I., and Beekman, F. J., “Review and current status of SPECT scatter correction,” Physics in Medicine and Biology 56(14), R85–R112 (2011).
  * [2] Garcia, E. V., “SPECT attenuation correction: an essential tool to realize nuclear cardiology’s manifest destiny,” Journal of Nuclear Cardiology 14(1), 16–24 (2007).
  * [3] Li, T. and Mok, G., “Mis-registration effect between quantitative SPECT and CT on activity estimation for targeted radionuclide therapy dosimetry,” Journal of Nuclear Medicine 59(supplement 1), 1746–1746 (2018).
  * [4] Sakoshi, M., Matsutomo, N., Yamamoto, T., and Sato, E., “Effect of misregistration between SPECT and CT images on attenuation correction for quantitative bone SPECT imaging,” Nihon Hoshasen Gijutsu Gakkai Zasshi 74(5), 452–458 (2018).
  * [5] Saleki, L., Ghafarian, P., Bitarafan-Rajabi, A., Yaghoobi, N., Fallahi, B., and Ay, M. R., “The influence of misregistration between CT and SPECT images on the accuracy of CT-based attenuation correction of cardiac SPECT/CT imaging: Phantom and clinical studies,” Iranian Journal of Nuclear Medicine 27(2), 63–72 (2019).
  * [6] Yoshikawa, K., Chono, T., Asanuma, O., and Yama, N., “Impact of misregistration between the myocardial perfusion images and CT attenuation correction map on the uptake with 17 segments polar map,” Nihon Hoshasen Gijutsu Gakkai Zasshi 75(9), 892–899 (2019).
  * [7] Tsuchiya, K., Takahashi, I., Kawaguchi, T., Yokoi, K., Morimoto, Y., Ishitsu, T., Suzuki, A., Ueno, Y., and Kobashi, K., “Basic performance and stability of a CdTe solid-state detector panel,” Annals of Nuclear Medicine 24(4), 301–311 (2010).
  * [8] Suzuki, A., Takeuchi, W., Ishitsu, T., Tsuchiya, K., Morimoto, Y., Ueno, Y., Kobashi, K., Kubo, N., Shiga, T., and Tamaki, N., “High-sensitivity brain SPECT system using cadmium telluride (CdTe) semiconductor detector and 4-pixel matched collimator,” Physics in Medicine & Biology 58(21), 7715 (2013).
  * [9] Pan, T.-S., King, M. A., De Vries, D. J., and Ljungberg, M., “Segmentation of the body and lungs from Compton scatter and photopeak window data in SPECT: a Monte-Carlo investigation,” IEEE Transactions on Medical Imaging 15(1), 13–24 (1996).
  * [10] Núñez, M., Prakash, V., Vila, R., Mut, F., Alonso, O., and Hutton, B. F., “Attenuation correction for lung SPECT: evidence of need and validation of an attenuation map derived from the emission data,” European Journal of Nuclear Medicine and Molecular Imaging 36(7), 1076–1089 (2009).
  * [11] Zaidi, H. and Hasegawa, B., “Determination of the attenuation map in emission tomography,” Journal of Nuclear Medicine 44(2), 291–315 (2003).
  * [12] Censor, Y., Gustafson, D. E., Lent, A., and Tuy, H., “A new approach to the emission computerized tomography problem: simultaneous calculation of attenuation and activity coefficients,” IEEE Transactions on Nuclear Science 26(2), 2775–2779 (1979).
  * [13] Krol, A., Bowsher, J. E., Manglos, S. H., Feiglin, D. H., Tornai, M. P., and Thomas, F. D., “An EM algorithm for estimating SPECT emission and transmission parameters from emission data only,” IEEE Transactions on Medical Imaging 20(3), 218–232 (2001).
  * [14] McMillan, A. B. and Bradshaw, T. J., “AI-based data corrections for attenuation and scatter in PET and SPECT,” arXiv preprint arXiv:2108.12412 (2021).
  * [15] Hagio, T., Poitrasson-Rivière, A., Moody, J. B., Renaud, J. M., Arida-Moody, L., Shah, R. V., Ficaro, E. P., and Murthy, V. L., ““Virtual” attenuation correction: improving stress myocardial perfusion SPECT imaging using deep learning,” European Journal of Nuclear Medicine and Molecular Imaging , 1–10 (2022).
  * [16] Yang, J., Shi, L., Wang, R., Miller, E. J., Sinusas, A. J., Liu, C., Gullberg, G. T., and Seo, Y., “Direct attenuation correction using deep learning for cardiac SPECT: a feasibility study,” Journal of Nuclear Medicine 62(11), 1645–1652 (2021).
  * [17] Chen, X., Zhou, B., Shi, L., Liu, H., Pang, Y., Wang, R., Miller, E. J., Sinusas, A. J., and Liu, C., “CT-free attenuation correction for dedicated cardiac SPECT using a 3D dual squeeze-and-excitation residual dense network,” Journal of Nuclear Cardiology , 1–16 (2021).
  * [18] Shi, L., Onofrey, J. A., Liu, H., Liu, Y.-H., and Liu, C., “Deep learning-based attenuation map generation for myocardial perfusion SPECT,” European Journal of Nuclear Medicine and Molecular Imaging 47(10), 2383–2395 (2020).
  * [19] Chen, X., Zhou, B., Xie, H., Shi, L., Liu, H., Holler, W., Lin, M., Liu, Y.-H., Miller, E. J., Sinusas, A. J., et al., “Direct and indirect strategies of deep-learning-based attenuation correction for general purpose and dedicated cardiac SPECT,” European Journal of Nuclear Medicine and Molecular Imaging , 1–15 (2022).
  * [20] Barrett, H. H., Myers, K. J., Hoeschen, C., Kupinski, M. A., and Little, M. P., “Task-based measures of image quality and their relation to radiation dose and patient risk,” Physics in Medicine & Biology 60(2), R1 (2015).
  * [21] Barrett, H. H., Yao, J., Rolland, J. P., and Myers, K. J., “Model observers for assessment of image quality.,” Proceedings of the National Academy of Sciences 90(21), 9758–9765 (1993).
  * [22] Jha, A. K., Myers, K. J., Obuchowski, N. A., Liu, Z., Rahman, M. A., Saboury, B., Rahmim, A., and Siegel, B. A., “Objective task-based evaluation of artificial intelligence-based medical imaging methods: Framework, strategies, and role of the physician,” PET Clinics 16(4), 493–511 (2021).  Artificial Intelligence and PET Imaging, Part I.
  * [23] Jha, A. K., Bradshaw, T. J., Buvat, I., Hatt, M., Prabhat, K., Liu, C., Obuchowski, N. F., Saboury, B., Slomka, P. J., Sunderland, J. J., et al., “Nuclear medicine and artificial intelligence: Best practices for evaluation (the RELAINCE guidelines),” Journal of Nuclear Medicine (2022).
  * [24] Yu, Z., Rahman, M. A., Schindler, T., Gropler, R., Laforest, R., Wahl, R., and Jha, A., “AI-based methods for nuclear-medicine imaging: Need for objective task-specific evaluation,” Journal of Nuclear Medicine 61(supplement 1), 575–575 (2020).
  * [25] Li, K., Zhou, W., Li, H., and Anastasio, M. A., “Assessing the impact of deep neural network-based image denoising on binary signal detection tasks,” IEEE Trans. Med. Imag. 40(9), 2295–2305 (2021).
  * [26] Yu, Z., Rahman, M. A., and Jha, A. K., “Investigating the limited performance of a deep-learning-based SPECT denoising approach: an observer-study-based characterization,” in [Medical Imaging 2022: Image Perception, Observer Performance, and Technology Assessment ], 12035, SPIE (2022).
  * [27] Rahman, M. A., Zhu, Y., Clarkson, E., Kupinski, M. A., Frey, E. C., and Jha, A. K., “Fisher information analysis of list-mode SPECT emission data for joint estimation of activity and attenuation distribution,” Inverse problems 36(8), 084002 (2020).
  * [28] Yu, Z., Rahman, M. A., Schindler, T., Laforest, R., and Jha, A. K., “A physics and learning-based transmission-less attenuation compensation method for spect,” in [Medical Imaging 2021: Physics of Medical Imaging ], 11595, 276–284, SPIE (2021).
  * [29] Hudson, H. M. and Larkin, R. S., “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Transactions on Medical Imaging 13(4), 601–609 (1994).
  * [30] Merlin, T., Stute, S., Benoit, D., Bert, J., Carlier, T., Comtat, C., Filipovic, M., Lamare, F., and Visvikis, D., “CASToR: a generic data organization and processing code framework for multi-modal and multi-dimensional tomographic reconstruction,” Physics in Medicine & Biology 63(18), 185005 (2018).
  * [31] Ronneberger, O., Fischer, P., and Brox, T., “U-net: Convolutional networks for biomedical image segmentation,” in [International Conference on Medical image computing and computer-assisted intervention ], 234–241, Springer (2015).
  * [32] Leung, K. H., Marashdeh, W., Wray, R., Ashrafinia, S., Pomper, M. G., Rahmim, A., and Jha, A. K., “A physics-guided modular deep-learning based automated framework for tumor segmentation in PET,” Physics in Medicine & Biology 65(24), 245032 (2020).
  * [33] Yousefirizi, F., Jha, A. K., Brosch-Lenz, J., Saboury, B., and Rahmim, A., “Toward high-throughput artificial intelligence-based segmentation in oncological PET imaging,” PET Clinics 16(4), 577–596 (2021).
  * [34] Schlemper, J., Oktay, O., Schaap, M., Heinrich, M., Kainz, B., Glocker, B., and Rueckert, D., “Attention gated networks: Learning to leverage salient regions in medical images,” Medical Image Analysis 53, 197–207 (2019).
  * [35] Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R., “Dropout: a simple way to prevent neural networks from overfitting,” The Journal of Machine Learning Research 15(1), 1929–1958 (2014).
  * [36] Brown, S., Bailey, D. L., Willowson, K., and Baldock, C., “Investigation of the relationship between linear attenuation coefficients and CT hounsfield units using radionuclides for SPECT,” Applied Radiation and Isotopes 66(9), 1206–1212 (2008).
  * [37] Narayanan, M. V., King, M. A., Leppo, J., Dahlbert, S., Pretorius, P. H., and Gifford, H. C., “Optimization of regularization of attenuation and scatter-corrected/sup 99m/Tc cardiac SPECT studies for defect detection using hybrid images,” IEEE Transactions on Nuclear Science 48(3), 785–789 (2001).
  * [38] Zhang, Y., Brady, M., and Smith, S., “Segmentation of brain MR images through a hidden Markov random field model and the expectation-maximization algorithm,” IEEE Transactions on Medical Imaging 20(1), 45–57 (2001).
  * [39] Glorot, X. and Bengio, Y., “Understanding the difficulty of training deep feedforward neural networks,” in [Proceedings of the thirteenth international conference on artificial intelligence and statistics ], 249–256, JMLR Workshop and Conference Proceedings (2010).
  * [40] Kingma, D. P. and Ba, J., “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980 (2014).
  * [41] Frey, E. C., Gilland, K. L., and Tsui, B. M., “Application of task-based measures of image quality to optimization and evaluation of three-dimensional reconstruction-based compensation methods in myocardial perfusion SPECT,” IEEE Transactions on Medical Imaging 21(9), 1040–1050 (2002).
  * [42] He, X., Frey, E. C., Links, J. M., Gilland, K. L., Segars, W. P., and Tsui, B. M., “A mathematical observer study for the evaluation and optimization of compensation methods for myocardial SPECT using a phantom population that realistically models patient variability,” IEEE Transactions on Nuclear Science 51(1), 218–224 (2004).
  * [43] Li, X., Jha, A. K., Ghaly, M., Elshahaby, F. E. A., Links, J. M., and Frey, E. C., “Use of sub-ensembles and multi-template observers to evaluate detection task performance for data that are not multivariate normal,” IEEE Transactions on Medical Imaging 36(4), 917–929 (2017).
  * [44] Wollenweber, S., Tsui, B., Lalush, D., Frey, E., LaCroix, K., and Gullberg, G., “Comparison of hotelling observer models and human observers in defect detection from myocardial SPECT imaging,” IEEE Transactions on Nuclear Science 46(6), 2098–2103 (1999).
  * [45] Sankaran, S., Frey, E. C., Gilland, K. L., and Tsui, B. M., “Optimum compensation method and filter cutoff frequency in myocardial SPECT: a human observer study,” Journal of Nuclear Medicine 43(3), 432–438 (2002).
  * [46] Barrett, H. H. and Myers, K. J., [Foundations of image science ], John Wiley & Sons (2013).
  * [47] Metz, C. E., “Roc methodology in radiologic imaging.,” Investigative radiology 21(9), 720–733 (1986).
  * [48] Metz, C. E., Herman, B. A., and Shen, J.-H., “Maximum likelihood estimation of receiver operating characteristic (ROC) curves from continuously-distributed data,” Statistics in Medicine 17(9), 1033–1053 (1998).
  * [49] Rahman, M. A. and Jha, A. K., “Task-based assessment of binned and list-mode SPECT systems,” in [Medical Imaging 2021: Image Perception, Observer Performance, and Technology Assessment ], 11599, 1159904, SPIE (2021).
  * [50] Rahman, M. A., Laforest, R., and Jha, A. K., “A list-mode OSEM-based attenuation and scatter compensation method for SPECT,” in [2020 IEEE 17th International Symposium on Biomedical Imaging (ISBI) ], 646–650, IEEE (2020).
  * [51] Clarkson, E. and Kupinski, M., “Quantifying the loss of information from binning list-mode data,” JOSA A 37(3), 450–457 (2020).
  * [52] Jha, A. K., Barrett, H. H., Frey, E. C., Clarkson, E., Caucci, L., and Kupinski, M. A., “Singular value decomposition for photon-processing nuclear imaging systems and applications for reconstruction and computing null functions,” Physics in Medicine & Biology 60(18), 7359 (2015).
  * [53] Jha, A. K. and Frey, E. C., “Estimating ROI activity concentration with photon-processing and photon-counting SPECT imaging systems,” in [Medical Imaging 2015: Physics of Medical Imaging ], 9412, 178–188, SPIE (2015).
  * [54] Van Rikxoort, E. M. and Van Ginneken, B., “Automated segmentation of pulmonary structures in thoracic computed tomography scans: a review,” Physics in Medicine & Biology 58(17), R187 (2013).
  * [55] Fu, Y., Lei, Y., Wang, T., Curran, W. J., Liu, T., and Yang, X., “A review of deep learning based methods for medical image multi-organ segmentation,” Physica Medica 85, 107–122 (2021).


# Dynamic Cluster Analysis to Detect and Track Novelty in Network Telescopes

Kai Huang  Politecnico di Torino
Turin, Italy
<EMAIL_ADDRESS>Luca Gioacchini  Politecnico di Torino
Turin, Italy
<EMAIL_ADDRESS>Marco Mellia  Politecnico di Torino
Turin, Italy
<EMAIL_ADDRESS>Luca Vassio  Politecnico di Torino
Turin, Italy
<EMAIL_ADDRESS>
###### Abstract

In the context of cybersecurity, tracking the activities of coordinated hosts
over time is a daunting task because both participants and their behaviours
evolve at a fast pace. We address this scenario by solving a dynamic novelty
discovery problem with the aim of both re-identifying patterns seen in the
past and highlighting new patterns. We focus on traffic collected by Network
Telescopes, a primary and noisy source for cybersecurity analysis. We propose
a 3-stage pipeline: (i) we learn compact representations (embeddings) of hosts
through their traffic in a self-supervised fashion; (ii) via clustering, we
distinguish groups of hosts performing similar activities; (iii) we track the
cluster temporal evolution to highlight novel patterns. We apply our
methodology to 20 days of telescope traffic during which we observe more than
8 thousand active hosts. Our results show that we efficiently identify 50-70
well-shaped clusters per day, 60-70% of which we associate with already
analysed cases, while we pinpoint 10-20 previously unseen clusters per day.
These correspond to activity changes and new incidents, of which we document
some. In short, our novelty discovery methodology enormously simplifies the
manual analysis the security analysts have to conduct to gain insights to
interpret novel coordinated activities.

###### Index Terms:

Network telescope, clustering, dynamic cluster analysis, novelty detection

##  1 Introduction

In the context of network monitoring and cybersecurity, the ever-evolving
landscape presents a perpetual challenge, with the daily emergence of new
threats and vulnerabilities [22, 42]. Take botnets which have always attracted
the attention of network administrators and security experts. Botnets are
groups of machines compromised by the same malware acting under the control of
a central server (Command-and-Control paradigm [35]). These compromised
machines become part of a large system that the controller can remotely
instruct to engage in coordinated malicious activities.

Botnets constantly evolve, adapting their attack patterns and network scanning
behaviours to the latest vulnerabilities. For instance, the FritzFrog botnet
[16] started exploiting the Log4Shell vulnerability [43] to improve its
spreading strategy [9]. Similarly, the widely known Mirai botnet [4, 5]
evolved to at least three new variants only in 2023 [32, 33]. This dynamic
nature requires a proactive and collaborative cybersecurity approach for
effective mitigation and eradication [34].

_How to represent and detect such coordinated activities_ among senders is a
crucial aspect of cyber defence. Solutions based on Artificial Intelligence
(AI) have proven instrumental in identifying these coordinated groups of
senders, offering a proactive approach to cybersecurity. Many works represent
senders through traditional feature engineering approaches and process the
resulting datasets through (sparse) Autoencoders [25, 18, 21], traditional
Convolutional Neural Networks [36, 3] or (temporal) Graph Neural Networks
[13]. Other works identify analogies between sequences of packets (or flows)
and words in text documents. Hence, they represent senders borrowing
techniques from the Natural Language Processing (NLP) field [37, 8, 14, 15].
The detection of senders performing similar actions is then faced by
addressing a clustering problem through unsupervised techniques in the
embedding space [18, 14, 15].

However, _distinguishing clusters and activity patterns that are already known
from those that are novel_ is still an open challenge. Consider the clustering
results obtained from analysing data referring to two periods of time, e.g.,
yesterday and today. Dynamic Cluster Analysis (DCA) aims at identifying which
of today’s clusters are novel and which can be mapped to yesterday’s activity.
Among many evolutionary clustering techniques for data stream [26, 23, 27],
MONIC [40] is one of the most adopted solutions. It defines a framework for
modelling and tracking cluster transitions over time. In short, it defines a
set of possible evolutions to describe cluster dynamics and to match clusters
across temporal snapshots.

Clustering techniques are commonly used to analyse network telescope traffic,
with embeddings [14], [21], leveraging feature vectors [19], [6], or based on
graphs [39]. However, very few works monitor the evolution of activities and
detect novelties. The closest to our work is [21], in which the authors
proposed an optimal mass transport problem to detect temporal changes in
clusters, whereas their approach focuses on unveiling large-scale incidences
like Mirai onset in 2016 but lacks tracking the evolution of individual
clusters and identifying newly emerged activities.

In this paper, we address these challenges by solving a dynamic novelty
discovery problem focusing on Network Telescopes, or Darknets. They are
subnets that host neither production services nor clients [12] and only
observe received unsolicited traffic. They represent a privileged source of
information for network security and monitoring activities [11, 20]. Here, we
analyse 20 days of telescope traffic, i.e., 20 independent snapshots of data.
We automatically detect groups of senders exhibiting coordination in the
traffic they send and study the evolution and dynamics over time of such
groups.

Figure 1: 3-stage pipeline for dynamic novelty discovery.

In short, we start the known 2-stage pipeline [24, 17, 15] in which, given a
batch of data, a self-supervised upstream task (Stage 1) generates compact
representations of input data in a latent space (i.e. embeddings) without the
need for ground truth or prior knowledge on the data; a clustering algorithm
(Stage 2) identifies groups of senders exhibiting similar patterns on top of
the generated embeddings. Here, we introduce a novel DCA stage (Stage 3) where
we match the current clusters with previously seen ones. This allows us to
investigate evolving dynamics, and pinpoint new previously unseen events.

Our analysis shows that the proposed approach greatly simplifies the analysis
of the humongous amount of data a network telescope collects. Each day, we
identify 50-70 clusters of coordinated senders. We immediately match 60-70% of
them to previously known clusters and identify only 10-20 new clusters, each
collecting hundreds of senders that exhibit a clear pattern. We manually
analyse some of them, identifying some periodic network scans, the sudden
emergence of attacks, morphing patterns in a group of senders, etc.

In summary, the key contributions of this work are as follows: We proposed a
DCA approach for detecting and tracking coordinated scanning activities on the
internet by presenting a novel framework based on MONIC. We apply our approach
to 20 days of real-world data collected from a /24 network telescope. The
results show our methodology successfully detects novel activities, limiting
human effort. We also manually analyse some of the identified cases and report
our findings.

All in all, our proposed pipeline demonstrates how modern machine-learning
approaches can support the analysis of traffic for cybersecurity goals. In the
future, we plan to move into the real-time implementation of the system and
continuous analysis.

##  2 Network Traffic Analysis Pipeline

Figure 1 shows the 3-stage pipeline we detail in the following. Starting from
the data collection, we summarise the self-supervised construction of
embedding (1) and the unsupervised cluster creation (2) that we borrow from
previous works. We then present the details of the new DCA process (3) we
introduce in this paper.

TABLE I: Traffic from active senders data characterisation. Total refers to distinct senders and ports. | 20 Days | Avg. Daily$\pm$std.
---|---|---
|  IPs |  Ports | IPs | Ports
TCP | 48,275 | 65,536 |  7,146$\pm$376 |  32,149$\pm$9,367
UDP | 5,825 | 3,897 |  1,413$\pm$105 |  413$\pm$105
ICMP | 930 | – |  151$\pm$42 | –
GRE | 1,396 | – |  286$\pm$154 | –
Total | 52,418 | 65,536 |  8,088$\pm$443 |  32,294$\pm$9,326

###  2.1 Network Telescope Sensor

In this work, we rely on data collected from the same /24 telescope as in [15]
situated in our campus network. We focus on 50 days of telescope traffic, from
2021-05-01 to 2021-06-19. In total, we observed more than 100 million packets
sent by 785 thousand senders. Here we consider senders sending more than 5
packets daily as “active”. After filtering out the inactive senders we are
left with about 130,000 active senders in the whole period. We use the first
30 days for bootstrapping the system and the remaining 20 days for testing it.
Figure 2 shows the total number of senders seen each day (green) (about
45,000). About 8,000 are _active_ each day (blue), with about 2,000 _new_
(red) daily senders appearing for the first time even after 20 days. In fact,
most senders are active for very few days as the Empirical Complementary
Cumulative Density Function (ECCDF) shows in Figure 3.

We detail some further statistics in Table I. In short, we face a very
variable scenario, with 52,418 active senders hitting the network telescope in
the last 20 days, of which each day $\approx$8,000 are active. They target all
TCP and thousands UDP ports, but also generate a significant fraction of ICMP
and GRE traffic.

###  2.2 Self-Supervised Upstream Task: Sender Embeddings

We generate sender embeddings through the state-of-the-art approach i-DarkVec
that we proposed [15], which relies on Word2Vec [29] and incremental training
to automatically create a vector representation for each sender sending
traffic to the traffic telescope in a given period. We assume that the reader
is familiar with Word2Vec and provide a brief overview of i-DarkVec. i-DarkVec
receives the sequence of packets as observed by the telescope in one day. Its
goal is to represent senders (identified by their IP addresses) so that
senders performing similar activities in the same period get projected in the
same portion of the space. For this, i-Darkvec extracts the sequences of
packets having the same protocol (ICMP/GRE/UDP/TCP) and set of destination
ports for TCP/UDP packets. Analogously to NLP, senders represent “words” and
their sequence represents “sentences”. We feed the generated sentences as
input to Word2Vec to map each sender to a vector in an $E$-dimension space.
i-DarkVec (Word2Vec) produces sender (word) embeddings such that senders
(words) co-occurring in sequences (sentences) appear close in the latent space
[15].

Figure 2: Daily senders on the 20-day dataset. Figure 3: ECCDF of number of
active days per sender on the 20-day dataset. Only 20% of the senders are
active for 3 days or more.

Formally, given a vocabulary of senders $V=\\{v_{1},v_{2},\ldots\\}$ and the
sequences of packets they send we map each entity ${v_{j}\in V\to
u_{j}\in\mathbb{R}^{E}}$ where $u_{j}$ is the embedding of $v_{j}$ in the $E$
dimensional space. The function $e:V\to\mathbb{R}^{E}$ is the embedding
function (i.e. Word2Vec) that we train by giving in input the sequences of
packets.

We consider each day as a new batch of data and we process it to extract the
embedding of each day.

###  2.3 Downtream Task: Unsupervised Clustering

We apply unsupervised downstream task on the produced embeddings using HDBSCAN
[28], a non-overlapping density-based clustering method. We employ the cosine
distance as the metric to measure the distance between sender embeddings.
Different from DBSCAN [10], HDBSCAN simplifies the parameter selection. It
only requires the minClusterSize, an intuitive and easy-to-choose parameter
that constrains the smallest size of clusters one wants.

In short, HDBSCAN effectively identifies clusters of points that lie in a
dense region of the space. It does not assume any shape for the cluster and
accommodates clusters with diverse internal densities. It extracts well-shaped
clusters, enhancing its robustness against noise. Points that fall in dense
regions with less than minClusterSize will be assigned the noise cluster.

Formally, at each discrete time window $t$, we partition the senders $V^{t}$
into clusters $X_{i}^{t}$, $V^{t}=\bigcup_{i}X_{i}^{t}$. We call the partition
$\zeta_{t}:=\\{X_{i}^{t}\\}_{i}$. The cluster $X_{-1}^{t}$ with the index
$i=-1$ is the noise cluster.

We assess the cluster quality by measuring the _Silhouette_ , $Sh\in[-1,1]$, a
measure of the cluster separation. Values close to 1 indicate perfect
separation of clusters. An average silhouette of over 0.7 is considered to be
strong, a value over 0.5 reasonable, and over 0.25 weak. Negative values are
an indication of bad clustering.

Given the embedding obtained for each batch of data $t$, we cluster active
senders to obtain $\zeta_{t}$.

(a) Number of detected clusters per day.

(b) ECDF of silhouette values.

(c) ECDF of cluster sizes.

Figure 4: Statistics of clustering results.

###  2.4 Dynamic Clusters Analysis

At the end of each time batch $t$ we obtain an independent clustering
$\zeta_{t}$. Hence, we need a way to track changes and detect novelties among
the clusters, i.e., to solve the DCA problem. MONIC [41] is a framework that
implements a DCA to model and track cluster transitions. It takes as an input
the clustering outcomes, e.g., $\zeta_{t}$ and $\zeta_{t+1}$ and defines a set
of transitions that identify possible changes involving multiple clusters.
Consequently, it offers valuable insights into the overall cluster dynamics.
To align the MONIC method with the IP embedding clusters scenario and to
enhance its suitability, we introduce some modifications. This adaptation
accounts for the presence of samples being clustered as noise and the
significant variation in samples for clustering (i.e., the active senders)
across different batches which are peculiar to our scenario.

Our DCA problem deals with the changes in members of clusters rather than
their spatial properties. Thus, the fundamental building block is the _set-
overlap_ operator defined as $OL(X,Y):=\frac{|X\cap Y|}{|X|}$, where $X$ and
$Y$ are two sets. Notice that the overlap operator is not symmetrical.

Given two clustering outcomes $\zeta_{t}$ and $\zeta_{t+1}$, we identify
various types of transitions according to the overlap among clusters. Let
$X_{i}$ be a cluster in $\zeta_{t}$ and $Y_{j}$ be a cluster in
$\zeta_{t+1}$.111For simplicity, we omit the index $t$ from the notation of
the clusters. Different from MONIC, we include the noise cluster in the
partitions, denoted as $X_{-1}$ and $Y_{-1}$.

Some senders active in $t$ may not be active in the following time slot, while
some senders would still active in $t+1$. The fraction of still present
senders is key to re-identify a cluster. Let us define the operator $A(X)$ as
the number of active elements of $X\in\zeta_{t}$ that are found in $t+1$:

$A(X):={\sum_{Y\in\zeta_{t+1}}OL(X,Y)}.$

Let $\tau_{0}\in[0.5,1]$ and $\tau_{1}\in(0,0.5)$ be two thresholds.

If $OL(X,Y)/A(X)\geq\tau_{0}$, we say that $X$ strongly matches $Y$. If
$OL(X,Y)/A(X)\geq\tau_{1}$, we say that $X$ loosely matches $Y$.222This is a
simplification w.r.t. MONIC original paper [41] in which there is an ageing
function.

First, we verify if a cluster is not anymore active in $t+1$:

  * •

Inactive: $A(X_{i})<\tau_{1}$; i.e., the fraction of active senders in $t+1$
is smaller than $\tau_{1}$.

Otherwise, we say the cluster $X_{i}$ is Active.

For active clusters, we define the following transformations:

  * •

Disappeared: $\forall
Y_{j}\in\zeta_{t+1}:\frac{OL(X_{i},Y_{j})}{A(X_{i})}<\tau_{1}$ or
$\frac{OL(X_{i},Y_{-1})}{A(X_{i})}\geq\tau_{0}$; i.e., there is no loosely
match with any other clusters; or there is a strong match with the noise
cluster.

  * •

Split: $\forall Y_{j}\in\zeta_{t+1}:\frac{OL(X_{i},Y_{j})}{A(X_{i})}<\tau_{0}$
and $\exists
Y_{j}\in\zeta_{t+1},j\neq-1:\frac{OL(X_{i},Y_{j})}{A(X_{i})}\geq\tau_{1}$;
i.e., there is not a strong match and there exists at least a loosely match in
$t+1$.

  * •

Survived: $\exists
Y_{j}\in\zeta_{t+1},j\neq-1:\frac{OL(X_{i},Y_{j})}{A(X_{i})}\geq\tau_{0}$ and
$\nexists k\neq i:\frac{OL(X_{k},Y_{j})}{A(X_{k})}\geq\tau_{0}$; i.e., $X_{i}$
is the only cluster that strongly matches another cluster
$Y_{j}\in\zeta_{t+1}$. Similarly, $X_{i}$ survives if $\exists k_{2}\neq
i:\frac{OL(X_{k_{2}},Y_{j})}{A(X_{k_{2}})}\geq\tau_{0}$ and
$OL(Y_{j},X_{i})\geq\tau_{0}$; i.e., there are other clusters that strongly
matches $Y_{j}$, but $X_{i}$ has a backward overlap larger than $\tau_{0}$
($X_{i}$ is the main component of $Y_{j}$).

  * •

Absorbed:
$Y_{j}\in\zeta_{t+1},j\neq-1:\frac{OL(X_{i},Y_{j})}{A(X_{i})}\geq\tau_{0}$ and
$\exists k\neq i:\frac{OL(X_{k},Y_{j})}{A(X_{k})}\geq\tau_{0}$; i.e., A
cluster $X_{i}$ strongly matches the cluster $Y_{j}$, but there is also
another cluster $X_{k}\in\zeta_{t}$ that strongly matches $Y_{j}$. in this
case, we require that the backward overlap of $X_{i}$ w.r.t. $Y_{j}$ is
limited, i.e. $OL(Y_{j},X_{i})<\tau_{0}$.

At last, we consider the case new clusters emerge. In details, we consider
$Y_{j}$ in $\zeta_{t+1}$ as:

  * •

Emerged: $\forall
X_{i}\in\zeta_{t},i\neq-1:\;\frac{OL(X_{i},Y_{j})}{A(X_{i})}<\tau_{0}$; i.e.,
there is no cluster in $t$ that strongly matches $Y_{j}$. Or $\forall
X_{i}\in\zeta_{t}:OL(Y_{j},X_{i})<\tau_{0}$; i.e., there are two or more
clusters in $t$ that strongly matches $Y_{j}$, but none of them is a major
component of $Y_{j}$.

The transitions Inactive and Emerged are not defined in the original
formulation of MONIC. In addition, we added the backward overlap to make sure
that if most of the points of $Y_{j}$ are present in $X_{i}$ then $X_{i}$
survived (otherwise it would have been absorbed).

In our dynamic cluster analysis, given a newly emerged cluster $Y$ at time
$t+1$, we perform a backward matching also with the past batches,
$b=[t-1,t-2,\ldots,0]$. For every past cluster $X\in\cup_{b}(\zeta_{b})$, we
check whether there exists some old $X$ which overlaps $Y$, i.e. if
$OL(X,Y)\geq\tau_{0}$. If there is not such a matching, we consider $Y$ as a
Novelty, i.e. we observe $Y$ for the first time in $t+1$.

TABLE II: Cluster evolution on the 20 days of test according to MONIC framework. | Labelled | Unknown
---|---|---
| Clusters† | Avg. IPs | Avg. Sh | Avg. purity | Clusters | Avg. IPs | Avg. Sh
Absorbed | 3.5% | 110 | 0.56 | 0.89 | 2.6% | 69 | 0.63
Split | 2.1% | 186 | 0.47 | 0.75 | 2.8% | 72 | 0.60
Disappeared | 2.0% | 19 | 0.70 | 0.93 | 3% | 16 | 0.69
Survived | 85.2% | 41 | 0.79 | 0.95 | 71.8% | 119 | 0.73
Inactive | 7.2% | 38 | 0.86 | 0.96 | 19.7% | 174 | 0.72
Total | 785 | 46 | 0.78 | 0.95 | 490 | 125 | 0.72
Emerged†† | 115 | 74 | 0.70 | 0.90 | 127 | 124 | 0.71
|  † Percentages exclude the clusters on the last day, since we do not know
their future transitions.
---
|  †† We detail characteristics of Emerged clusters (already included in the
Total row).
---

##  3 Ground Truth for Testing

To verify the performance of both the clustering and the dynamic analysis we
leverage a partial ground truth (GT) representing 13 classes of senders whose
coordination is known a priori. These are mostly set of IP addresses belonging
to well-known projects or botnets. In detail, we label Mirai-like senders by
checking the Mirai packet signature [7], and we label acknowledged internet
scanners publicly advertised in online
repositories333https://gitlab.com/mcollins_at_isi/acknowledged_scanners. We
mark all the remaining senders as _Unknown_ , ending up with a highly
unbalanced GT – e.g. thousands of senders exhibit the Mirai-like label, while
only hundreds or dozens belong to classes of acknowledged scanners.

To assign a cluster label, we employ a majority voting rule, i.e. we assign a
cluster the most frequently occurring label within samples in that cluster. We
also consider the purity of a cluster[38], i.e. the proportion of samples
within a cluster bearing the predominant label.

In general, we expect that senders belonging to the same GT class shall
consistently belong to the same cluster which survives over time. We
eventually expect subgroups of senders in a given class to perform different
scan activities444We already observed subsets of senders of a known class that
perform different scans [15].. We instead expect new clusters to Emerge,
Disappear or be Absorbed for senders in the unknown class.

##  4 Experimental Results

The size of the embeddings $E$ is chosen as 200, as in [15]. We set HDBSCAN
with minClusterSize = 10 and $\tau_{0}=0.65$ and $\tau_{1}=0.3$ for the DCA
algorithm. While there is room for future tuning of these parameters, we
empirically found these choices have little impact on the overall results.

###  4.1 Generic Clustering Results

As previously said, we consider batch size equal to one day and run HDBSCAN
for the 20 testing days in our dataset. For each day, we apply clustering only
for active senders.

Figure 4(a) shows the number of clusters found in each day. We observe a quite
variable setup, with from 53 to 69 clusters found each day. Considering the
quality of clustering, in Figure 4(b) we show the empirical cumulative
distribution function (ECDF) of the silhouette for the clusters identified in
the initial, mid, and final days of the dataset. We observe some variations of
cluster quality across different dates, but typically more than 80% of
clusters have $sh\geq 0.6$, i.e., most of the clusters are well-separated and
compact. To complete the clustering characterization, Figure 4(c) shows the
ECDF of cluster sizes. Most of them have less than 100 senders but some tops
to 1,000 samples (notice the log x-scale).

In short, the clustering step produces few and well-shaped clusters, reducing
by two order of magnitude the entities to analyse. The question remains which
of the cluster present at batch $t+1$ can be linked to clusters seen in the
past batches. We answer to this question next.

(a) Cluster transitions per day.

(b) Emerged clusters per day.

Figure 5: Dynamic cluster analysis. Figure 6: Evolution of all activities.
Lines link clusters that split (dotted green line) or get absorbed (orange
dotted line).

(a) Shadowserver clusters.

(b) Censys clusters.

Figure 7: Examples of evolution of some clusters of known classes with DCA
framework.

###  4.2 Overall Results of DCA

In Figure 5, we evaluate cluster evolution with the outcomes of the DCA. Every
day, after obtaining new clustering, we backtrack changes in activities from
the previous days. Figure 5(a) shows the breakdown of clusters among Survived,
Inactive, Absorbed, Disappeared and Split. We observe that the majority of
clusters (40 to 50) survive, indicating most of the previously identified
coordinated activities persist. There is no need to analyse these new clusters
if one did it in the past.

A small portion of clusters becomes inactive (blue), i.e., their senders are
not present in $t+1$ anymore. Very few clusters occurred in some other
transition, most of which correspond to senders being Absorbed by larger
clusters (green) or getting Split (cyan), i.e., their behaviour is not part of
a strong coordinated pattern. By manually inspecting some of these cases, we
observe mostly noisy senders which send few packets but do not follow a
particular pattern.

Each day, the DCA algorithm identifies also some Emerged clusters (7-20 each
day). Figure 5(b) details this reporting those that will Survive in the next
days and those that will never be re-identified in the future. About half of
them would survive for at least one more day. The analysts should give
priority to checking these potentially new (and persistent) clusters that
could represent new anomalies and threats.

Using the DCA we identify 216 different coordinated activities, we can label
about 37.5% of them. Figure 6 provides an overview of the overall pattern of
these. On the x-axis, we have the time. The y-axis reports the cluster ID (in
increasing order). When a new cluster Emerges for the first time, we assign it
a new ID. If a cluster Survives, we report it by assigning the previous ID and
reporting it with a dash. When a cluster disappears, we mark the last time we
saw it by a corresponding ending event (Absorbed or Disappeared). In the case
of Split, we link the new clusters to the originating one (orange lines). Same
for Absorbed (green lines).

Some consideration holds: 1) some clusters are identified on the first day and
survive throughout the 20 days (bottom part of the plot); 2) some emerge later
and then survive for multiple days (linear growth followed by horizontal
dashes); 3) some become inactive for several days, then re-appear
(intermittent dash patterns). Only a few clusters Split or are Absorbed.

For completeness, we detail the transitions with labelled and unknown clusters
in Table II. As expected, clusters with labels are more stable because they
involve senders controlled by research projects or security crawlers that
exhibit regular and constant activities. The DCA can easily match most of them
as Surviving (85.2%). Notice 115 clusters Emerge555We do not consider the
clusters at the first day as Emerged., a small fraction compared to the 785
that we observe in the whole 20 days. Notice that these labelled clusters
exhibit high purity and silhouette, which confirms accurate clustering. A
small portion of clusters have low purity and tend to Split more easily. For
instance, some noise points are initially clustered together labelled points.
In the next batch, the clustering step may remove those noise points (which
may generate a separate cluster, resulting in a Split).

Focus now on unknown clusters – right part of the Table II. They manifest
greater dynamism with nearly 20% of them becoming inactive (w.r.t. 7.2% for
labelled clusters). Only 72% of clusters survive for more than one day (w.r.t.
85%). We also witness more Emerging clusters (127 Emerged over a total of 490
clusters identified over 20 days). Interestingly, unknown clusters are
generally larger than labelled ones. This can be linked to possible botnet
activities which can involve a larger number of senders than research projects
(which rely on limited and fixed senders acting as scanners). In fact, the
diversity of IP addresses and subnets found in these clusters confirms this
assumption.

###  4.3 Example of Labelled Cluster Evolution

To give the intuition of the goodness of the DCA process, we focus on the
clusters of senders belonging to two well-known acknowledged security engines:
Shadowserver [2] and Censys [1]. We know these projects keep periodically
scanning the entire IPv4 address space looking for vulnerable services. They
generate a continuous and regular stream of packets and we expect them to be
consistently present and well-identified by i-Darkvec and the clustering step,
while we expect the DCA to keep strongly matching current clusters with past
ones (Survived).

Focus on Shadowserver first. In Figure 7(a) we show the evolution of labelled
clusters with the majority of senders being known Shadowserver scanners. We
identify 16 clusters on the first day (we report them using their unique
cluster ID as it appears in Figure 6). This reflects the fact that different
groups of Shadowserver’s senders perform different scan activities. In the
following days, the DCA analysis shows a very stable picture which reflects
the periodic scan Shadowserver’s senders do – i.e. all of the clusters Survive
for 20 days. There are a few exceptions: clusters 9 and 10 which are absorbed
together into cluster 16 for some days. By manually verifying these cases, we
notice in fact that senders in these two clusters generate less traffic on
those days, and our pipeline reflects this.

(a) Online services scanning cluster.

(b) GRE backscattering cluster.

(c) 9530/TCP-scan clusters evolution.

Figure 8: Patterns of activity for senders of some example clusters.

Focus now on Censys in Figure 7(b). We identify 20 clusters, all labelled as
Censys. Here, some clusters become Inactive (missing dash) and reappear from
time to time. Some are Absorbed into clusters with unlabelled IPs (red dashed
lines in the plot), Split or even Disappear into noise. By manually checking
some of these events, we confirm that Censys scanners are split into groups
performing different scan activities. Their activity is less constant in time.
In fact, some senders sometimes disappear for days or portions of a day. This
alters the sequence of packets i-Darkvec sees, modifying the host embedding
and making these more similar to (and thus clustered with) other intermittent
or noisy senders. The DCA pinpoints these changes. When those senders become
active again, i-Darkvec and the clustering assign them to their original
cluster and the DCA matches it with its last occurrence.

These two examples show the robustness of our DCA methodology in successfully
clustering and matching senders performing similar activities.

##  5 Novelty Examples

We manually analyse some of the clusters (mostly unknown) and look for
insights about their activities. We leverage external information offered by
DNS and whois and summarise some interesting findings.

###  5.1 Detailed Examples

SIP and RDP scanners in Google Cloud Figure 8(a) shows the patterns of senders
in Cluster 27 that DCA identifies in the first day. The $x$-axis represents
the time, and the $y$-axis the senders, ordered by time of appearance; a point
marks the time in which the given sender sends a packet. The bottom 90 senders
are from three /24 subnets located in different Google Cloud IP address
ranges. Their activities show evident temporal patterns (see the bottom of the
figure), which persist for 20 days. These senders target port 5060/TCP (SIP
initialization protocol) and port 3389/TCP (used for Windows Remote Desktop -
RDP). Some of them are reported in some database for unauthorized connection
attempts. This evidence strongly suggests some malicious scan activity from
actors leveraging the Google Cloud. Senders 90-150 in the top part of the
Figure are less regular but perform scans to the same ports with a similar
pattern. Thus the clustering algorithm correctly places them in the same
cluster (C27). Their IP are like bots scattered in random networks and also
reported for abuse in some database. These senders could be controlled by
other entities, or be part of a distributed botnet the malicious actor
controls as well.

GRE DDoS attack Figure 8(b) shows the activity pattern exhibited by Cluster 62
(C62) senders over time. The $y$-axis shows the senders of this cluster in
order of appearance. Here we are also reporting ($x$-axis) what happens during
some prior days of bootstrapping. This cluster suddenly emerges and survives
for 20 days. In the first days, we observe more than 700 senders. The number
keeps shrinking, with about 100 senders still present on the last day. Around
30% of the senders belong to the domain name rr.com, a webmail provider. Most
of the traffic is Generic Routing Encapsulation (GRE) protocol. Senders hit
our network telescope on 2021-05-30, right before our test period and become
very active in the following few days. Manually checking the packet trace, we
conclude these are back-scattering packets sent by victims of a Distributed
Denial of Service (DDoS) attack where rr.com senders are the victims. The
attacker sends GRE packets to rr.com servers with randomly spoofed IP sender
addresses. By chance, these happen to be one of our /24 darknet addresses. The
victim then replies with a GRE error message, sending packets to our darknet
IP addresses. The decrease in senders is likely due to some mitigation action
the victim put in place.

9530/TCP scan clusters Figure 8(c) shows the evolution of Clusters 105, 176
and 187. These senders mostly target port 9530/TCP, which is reported as a
vulnerability (backdoor) of some surveillance devices [31]. We identify C105
firstly on 2021-06-03 (first day shown on the $x$-axis) with only 15 senders
(11 of them are Mirai-like senders). This cluster expands significantly up to
about 800 senders on 06-12 when suddenly we observe a lot of packets (see the
dense red pattern). Most of those senders go suddenly silent the day after,
and only a small portion remains still active. The DCA fails to match the new
cluster C176 with C105 (because $OL(C_{105},C_{176})/A(C_{105})<\tau_{0}$) as
most senders of C105 are inactive now). We thus identify C176 as newly emerged
and classify it as a novelty. Most of the inactive senders suddenly become
active again in 06-14 and are now matched with C176 (see the dense green
pattern). On 06-15, they suddenly become inactive, and the DCA identifies the
few remaining active as a new cluster – C187, whose members grow in the next
few days.

This example shows our approach successfully identifies the dynamic. Yet, the
extreme on/off/on pattern challenges our ability to match newly emerged
clusters with previous ones. This shows the complexity of possible scenarios,
and opens the opportunity to improve the DCA.

###  5.2 Other Examples

Mongolian /24 scanner (C0) The cluster appears for the first time on
2021-05-31 then became inactive for several days until the 2021-06-13 when it
survives for three days. In total, we observe 256 senders in this cluster, 253
of which belong to the same subnet 180.149.126.0/24 located in Mongolia. These
senders perform TCP scans to some very specific ports for the whole 48 hours
of the weekend in local time. Some of these addresses have been reported in
the past as a possible network
scanner666https://www.abuseipdb.com/check/180.149.126.40.

Targeted port scan (C7) The cluster survives throughout the 20 days. The
cluster is composed of 21 senders distributed in random networks (only four of
them are from the same /24 subnet). All senders were performing scans at a
high frequency for 20 days. However, on the evening of 06-07 all scan
activities paused for about 12 hours, which reflects an evident coordination.
Senders constantly target the same 35 ports, all well-known for vulnerable
services, e.g. port 1900 used by SSDP[30]. This clearly hints at a malicious
coordinated activity. Even if only 21 senders are involved, our methodology
could pinpoint and highlight it.

New Shadowserver senders (C20) The cluster survives for the whole 20 days,
being remarkably stable. It includes about 200 senders from the same /24
network 64.62.197.0/24. A reverse DNS lookup reports these belong to
Shadowserver, and their temporal patterns are consistent with the other known
Shadowserver scanners. Hence, we confirm a new group of Shadowserver senders
was added to their infrastructure, but not included in the publicly available
GT.

Botnet targeting port 23 (C61) The cluster is quite large and dynamic. Its
size varies dramatically day by day, from around 500 to more than 1,100. Each
day, about half of the senders in that cluster became inactive and hundreds of
new senders joined. Nevertheless, the cluster is always identified as
surviving. In total, we observe nearly 7 thousand senders in the cluster, with
more than 3,500 IPs just appearing once. They target mostly the 23/TCP port.
This suggests these senders are part of a large botnet. Interestingly, on the
first day, 90% of senders are labelled Mirai-like. The cluster purity
decreases from 0.9 to about 0.5 so that since day 13 the majority of senders
are unknown. After inspecting the traffic of these new senders, we find that
most of them exhibit the Mirai fingerprint. These are new senders that were
not present in the dataset we used to build the GT. This testifies the Mirai-
like botnet keeps growing, and that our approach can simplify the tracking of
new bots appearing.

##  6 Conclusions

In this paper, we presented a preliminary exploration to track the evolution
of coordinated activities and detect novelty from traffic observed via network
telescopes. We proposed a dynamic clusters analysis approach for activities in
network telescopes inspired by the MONIC framework.

Experimental results show that our methodology can: (i) reduce human effort to
analyse newly observed groups of coordinates senders; (ii) provide deeper
insights into activities, highlight changes and incidence by tracking the
clusters evolution.

Future development includes improving the robustness of the approach by
defining improved strategies to match clusters at distant times and thoroughly
exploring the sensitivity to parameters. A promising extension of the proposed
methodology is to leverage other features in addition to those extracted by
i-DarkVec, introducing criteria on priority to check new clusters and
automatically highlight changes within a cluster. Furthermore, we are looking
forward to deploying the system in real-time and performing continuous
analysis.

## Ethical Considerations

Our work focuses solely on technical advancements in utilising dynamic cluster
analysis to enhance cybersecurity and deepen the knowledge of traffic observed
in network telescopes. We only employ passive measurement techniques, this
means that we do not engage with or influence any entities we measure. No
further ethical considerations apply to our research, as our primary objective
is to improve the effectiveness of network security measures through
innovative AI solutions.

## Acknowledgements

This work has been partly funded by the project SERICS (SEcurity and RIghts In
the CyberSpace - PE00000014) under the MUR National Recovery and Resilience
Plan funded by the European Union, and the xInternet (eXplainable Internet -
20225CETN9) projects - funded by European Union - Next Generation EU within
the PRIN 2022 program (D.D. 104 - 02/02/2022 Ministero dell’Università e della
Ricerca). This manuscript reflects only the authors’ views and opinions and
the Ministry cannot be considered responsible for them.

## References

  * [1] Censys, 2021.  https://censys.io/.
  * [2] The shadowserver foundation, 2021.  https://www.shadowserver.org/.
  * [3] G. Aceto, D. Ciuonzo, A. Montieri, and A. Pescapè.  Mimetic: Mobile encrypted traffic classification using multimodal deep learning.  Computer Networks, 2019.
  * [4] Antonia Affinito, Stefania Zinno, Giovanni Stanco, Alessio Botta, and Giorgio Ventre.  The evolution of mirai botnet scans over a six-year period.  Journal of Information Security and Applications, 2023.
  * [5] Manos Antonakakis, Tim April, Michael Bailey, Matt Bernhard, Elie Bursztein, Jaime Cochran, Zakir Durumeric, J Alex Halderman, Luca Invernizzi, Michalis Kallitsis, et al.  Understanding the mirai botnet.  In 26th USENIX security symposium (USENIX Security 17), 2017.
  * [6] Olivier Cabana, Amr M Youssef, Mourad Debbabi, Bernard Lebel, Marthe Kassouf, Ribal Atallah, and Basile L Agba.  Threat intelligence generation using network telescope data for industrial control systems.  IEEE Transactions on Information Forensics and Security, 2021.
  * [7] J. Ceron, K. Steding-Jessen, C. Hoepers, L. Granville, and C. Margi.  Improving iot botnet investigation using an adaptive network layer.  Sensors, 2019.
  * [8] D. Cohen, Y. Mirsky, Y. Elovici, R. Puzis, M. Kamp, T. Martin, and A. Shabtai.  DANTE: A framework for mining and monitoring darknet traffic.  2020\.
  * [9] Ori David.  Frog4Shell — FritzFrog Botnet Adds One-Days to Its Arsenal, 2024.  Akamai.
  * [10] Martin Ester, Hans-Peter Kriegel, Jörg Sander, Xiaowei Xu, et al.  A density-based algorithm for discovering clusters in large spatial databases with noise.  In kdd, 1996.
  * [11] C. Fachkha, E. Bou-Harb, and M. Debbabi.  Inferring Distributed Reflection Denial of Service Attacks from Darknet.  Computer Communications, 2015.
  * [12] C. Fachkha and M. Debbabi.  Darknet as a Source of Cyber Intelligence: Survey, Taxonomy, and Characterization.  Commun. Surveys Tuts., 2016.
  * [13] Luca Gioacchini, Andrea Cavallo, Marco Mellia, and Luca Vassio.  Exploring temporal gnn embeddings for darknet traffic analysis.  In Proceedings of the 2nd on Graph Neural Networking Workshop 2023, 2023.
  * [14] Luca Gioacchini, Luca Vassio, Marco Mellia, Idilio Drago, Zied Ben Houidi, and Dario Rossi.  DarkVec: Automatic analysis of darknet traffic with word embeddings.  In Proceedings of the 17th International Conference on emerging Networking EXperiments and Technologies, 2021.
  * [15] Luca Gioacchini, Luca Vassio, Marco Mellia, Idilio Drago, Zied Ben Houidi, and Dario Rossi.  i-DarkVec: Incremental Embeddings for Darknet Traffic Analysis.  ACM Transactions on Internet Technology, 2023.
  * [16] Ophir Harpaz.  FritzFrog: A New Generation of Peer-to-Peer Botnets, 2022.  Akamai.
  * [17] Z.B. Houidi, R. Azorin, M. Gallo, A. Finamore, and D. Rossi.  Towards a systematic multi-modal representation learning for network data.  In Proceedings of the 21st ACM Workshop on Hot Topics in Networks, 2022.
  * [18] J. Höchst, L. Baumgärtner, M. Hollick, and B. Freisleben.  Unsupervised traffic flow classification using a neural autoencoder.  In 2017 IEEE 42nd Conference on Local Computer Networks (LCN), 2017.
  * [19] Félix Iglesias and Tanja Zseby.  Pattern discovery in internet background radiation.  IEEE Transactions on Big Data, 5(4):467–480, 2017.
  * [20] M. Jonker, A. King, J. Krupp, C. Rossow, A. Sperotto, and A. Dainotti.  Millions of Targets UnderAttack: A Macroscopic Characterization of the DoS Ecosystem.  In Proceedings of the ACM SIGCOMM Internet Measurement Conference, 2017.
  * [21] M. Kallitsis, R. Prajapati, V. Honavar, D. Wu, and J. Yen.  Detecting and interpreting changes in scanning behavior in large network telescopes.  IEEE Transactions on Information Forensics and Security, 2022.
  * [22] Pooja Kumari and Ankit Kumar Jain.  A comprehensive study of DDoS attacks over IoT network and their countermeasures.  Computers & Security, 2023.
  * [23] Hui Liu, Aihua Wu, Mingkang Wei, and Chin-Chen Chang.  Skdstream: a dynamic clustering algorithm on time-decaying data stream.  EURASIP Journal on Wireless Communications and Networking, 2022.
  * [24] Y. Liu, M. Jin, S. Pan, C. Zhou, Y. Zheng, F. Xia, and S.Y. Philip.  Graph self-supervised learning: A survey.  IEEE Transactions on Knowledge and Data Engineering, 2022.
  * [25] M. Lotfollahi, R.S.H. Zade, M.J. Siavoshani, and M. Saberian.  Deep packet: A novel approach for encrypted traffic classification using deep learning, 2017.
  * [26] Mahmoud A Mahdi, Khalid M Hosny, and Ibrahim Elhenawy.  Scalable clustering algorithms for big data: A review.  IEEE Access, 2021.
  * [27] José Maia, Carlos Alberto Severiano Junior, Frederico Gadelha Guimarães, Cristiano Leite de Castro, André Paim Lemos, Juan Camilo Fonseca Galindo, and Miri Weiss Cohen.  Evolving clustering algorithm based on mixture of typicalities for stream data mining.  Future Generation Computer Systems, 2020.
  * [28] Leland McInnes, John Healy, and Steve Astels.  HDBSCAN: Hierarchical Density Based Clustering.  The Journal of Open Source Software, 2017.
  * [29] T. Mikolov, K. Chen, G. Corrado, and J. Dean.  Efficient estimation of word representations in vector space.  arXiv, 2013.
  * [30] n.a.  Stupidly simple ddos protocol (ssdp) generates 100 gbps ddos, 2017.  Cloudflare,https://blog.cloudflare.com/ssdp-100gbps/.
  * [31] n.a.  Full disclosure: 0day vulnerability (backdoor) in firmware for Xiaongmai-based DVRs, NVRs and IP cameras, 2020.  Habr.
  * [32] n.a.  Mirai Botnet’s New Wave: hailBot,kiraiBot, catDDoS, and Their Fierce Onslaught, 2023.  NSFOCUS.
  * [33] n.a.  New Mirai-based Botnet Variants Emerge, 2023.  SOCRadar.
  * [34] n.a.  Qakbot Malware Disrupted in International Cyber Takedown, 2023.  Office of Public Affairs – U.S. Department of Justice.
  * [35] Sneha Padhiar, Aayushyamaan Shah, and Ritesh Patel.  The Hidden Enemy: A Botnet Taxonomy.  In Proceedings of the International Conference on Paradigms of Computing, Communication and Data Sciences: PCCDS 2022, 2023.
  * [36] S. Rezaei and X. Liu.  How to achieve high classification accuracy with just a few labels: A semi-supervised approach using sampled packets, 2018.
  * [37] M. Ring, A. Dallmann, D. Landes, and A. Hotho.  IP2Vec: Learning Similarities Between IP Addresses.  In 2017 IEEE International Conference on Data Mining Workshops (ICDMW), 2017.
  * [38] Hinrich Schütze, Christopher D Manning, and Prabhakar Raghavan.  Introduction to information retrieval, volume 39.  Cambridge University Press Cambridge, 2008.
  * [39] Francesca Soro, Mauro Allegretta, Marco Mellia, Idilio Drago, and Leandro M Bertholdo.  Sensing the noise: Uncovering communities in darknet traffic.  In 2020 Mediterranean Communication and Computer Networking Conference (MedComNet), pages 1–8. IEEE, 2020.
  * [40] Myra Spiliopoulou, Eirini Ntoutsi, Yannis Theodoridis, and Rene Schult.  MONIC and Followups on Modeling and Monitoring Cluster Transitions.  In Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2013, Prague, Czech Republic, September 23-27, 2013, Proceedings, Part III 13, 2013.
  * [41] Myra Spiliopoulou, Irene Ntoutsi, Yannis Theodoridis, and Rene Schult.  Monic: modeling and monitoring cluster transitions.  In Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 706–711, 2006.
  * [42] Nan Sun, Ming Ding, Jiaojiao Jiang, Weikang Xu, Xiaoxing Mo, Yonghang Tai, and Jun Zhang.  Cyber threat intelligence mining for proactive cybersecurity defense: A survey and new perspectives.  IEEE Communications Surveys & Tutorials, 2023.
  * [43] Yue Wang, Md Abul Bashar, Mahinthan Chandramohan, and Richi Nayak.  Exploring Topic Models to Discern Cyber Threats on Twitter: A Case Study on Log4Shell.  Intelligent Systems with Applications, 2023.