Multiple myeloma (MM) is a heterogeneous disease regarding symptoms, tumor genetics, and outcome.1 Precise diagnostics to classify disease subtypes and to predict the individual course of disease is of utmost importance to guide treatment accordingly, a principle known as precision oncology. In addition to serological and urine parameters, imaging and bone marrow (BM) biopsy results play decisive roles in diagnosis of monoclonal plasma cell disorders. Parameters received from unguided BM biopsy at the posterior iliac crest as plasma cell infiltration (PCI)2–7 or cytogenetic aberrations8–11 have also proven value as biomarkers and consequently are now used for staging, risk stratification, and response assessment.9,10,12–15 However, tumor load distribution16,17 and genomic aberrations18–20 can be spatially heterogeneous. Invasive biopsies have the disadvantage that they cannot be performed both multifocally and in frequent repetition, which would be necessary to precisely evaluate and monitor the tumor load and capture the complete genomic landscape in each patient. Here, whole-body magnetic resonance imaging (wb-MRI) is advantageous, as it allows to investigate the complete BM of a patient noninvasively and to capture information on spatial distribution and local characteristics of tumor manifestations. Consequently, a method to predict local biopsy results from MRI without multiple invasive biopsies would be of great value, if a link between local imaging findings and local histology or local genetic findings could be established.

Radiomics is a new image analysis approach, characterizing a structure from imaging by calculating hundreds of mathematically defined radiomics features to quantify signal intensity, shape, and texture of the structure.21,22 For various oncologic entities, it has been reported that radiomics can predict tumor tissue characteristics, as histologic or genetic results.23 For patients with MM, it was recently demonstrated that the combination of an nnU-Net24 for segmentation and a subsequent radiomics analysis21,22 allows to automatically obtain an objective, in-depth characterization of the whole BM automatically from wb-MRI.25

Given this clinical, biological, and technical background, the purpose of this study was to develop and test an automated deep learning and radiomics image analysis framework, which analyzes the pelvic BM from whole-body MRIs to predict the local tumor tissue characteristics PCI and cytogenetic aberrations from routine BM biopsy at the iliac crest.

METHODS Study Design and Algorithmic Concept

This study was designed as a retrospective study using multicentric data sets to establish and test a methodological concept to perform automatic tissue analysis from routine clinical MRIs and predict BM biopsy results in patients with monoclonal plasma cell disorders. To achieve this goal, an automated multistep pipeline was established, which included automated pelvic BM segmentation by a deep learning algorithm,24 image normalization26 and resampling, radiomics feature calculation,27 and parameter prediction based on a machine learning model (Fig. 1). This study was approved by the institutional review board Medical Faculty, University of Heidelberg (Germany); S-537/2020 with waiver of informed consent. The acquisition of imaging and clinical data was performed between 2008 and 2021. Specific planning, data annotation, establishing algorithms, and data analysis for this project were performed between 2019 and 2022.

FIGURE 1: Multistep workflow for automated image analysis and prediction of biopsy results. Step 0: A coronal T1-weighted turbo spin echo sequence is used as input for the processing pipeline. Step 1: The nnU-Net, a deep learning-based segmentation algorithm,24 performs automatic segmentation of the right hip bone and the left hip bone and marks each side individually. In addition, the medial parts of the piriformis muscles are also automatically segmented (both sides receive the same label). Step 2: Images are normalized to the mean signal intensity of the piriformis muscle,26 and resampled to a uniform geometry. Step 3: 260 radiomics features (91 first order and 169 texture features) are calculated by the MITK phenotyping27 radiomics toolbox. Step 4: Machine-learning models predict the target parameters bone marrow plasma cell infiltration, cytogenetic risk status, and cytogenetic aberrations from the radiomics features. T1w TSE, T1-weighted turbo spin echo sequence.Patient Cohorts and Data Sets

There was no prespecified sample size for this study. As recommended, we aimed to assemble an overall data set with a large number of samples so that machine learning methods can be reasonably applied, with a balanced distribution of patients with low and high PCI, and which allows for independent, external, multicentric testing of the final algorithms. Data from center 1 were derived from 5 different subsets, which had in part been included in other studies. Data from center 1 were used for training of the algorithms, and a hold-out subset was reserved for independent, internal testing (internal test set). Data from center 2 were used for external testing: once in a subset of MRIs with a homogenized MRI protocol and high imaging quality (center 2, high-quality test set), and once in a subset that did not fulfill these criteria (center 2, other test set). Finally, a multicentric data set was used to evaluate the performance in a very heterogenous data set (multicenter test set). Details on the inclusion and exclusion process including the respective flow-charts are reported in Supplemental Digital Content 1, https://links.lww.com/RLI/A806. An overview over all data sets is shown in Figure 2.

FIGURE 2:

Overview of study cohorts and data sets. Data from center 1 were used to train the automatic segmentation algorithm and the radiomics algorithms, and an independent subset of 59 patients was preserved for internal testing of the resulting pipeline. aPatients without concomitant bone marrow biopsy were assigned to the training set for the nnU-Net only. Patients with concomitant bone marrow biopsy were split by date into training and test sets as described in detail in the supplements. bThis training part was then used both for training of the nnU-Net and the radiomics algorithms. Then, the resulting pipeline was tested on the internal test set. cBefore external testing, the radiomics algorithms were retrained on all internal data. Then, the radiomics algorithms were tested on 3 external test sets: (1) the high-quality test set from center 2 that had been acquired with a harmonized MRI protocol at the newest MRI scanner and was free of imaging artifacts; (2) all data from center 2 that did not meet these criteria (center 2, other test set); and (3) the multicenter test set that included data from a wide variety of scanners from 6 different centers.

Imaging

Coronal T1-weighted (T1w) turbo spin echo images that had been acquired with different MRI scanners and different MRI sequence parameters were included. Detailed information on all scanners and sequence parameters is reported in Supplemental Digital Content 2, https://links.lww.com/RLI/A807.

Image Segmentation and Training and Testing of nnU-Nets

The BM of the right and left hip bone and the medial part of the piriformis muscle were segmented on coronal T1w images. An initial subset of 127 MRIs was annotated manually to train a first nnU-Net,24 which was then used to automatically presegment the rest of the training data. After manual improvement by 1 experienced rater (7 years of experience in MRI segmentation), the final nnU-Net was trained on all training data (470 cases). Testing of the nnU-Net was performed in 3 independent test sets with overall 37 manually segmented cases as reference. Further details are reported in Supplemental Digital Content 3, https://links.lww.com/RLI/A808.

Radiomics Analysis

Before feature extraction, all images were resampled to a uniform voxel spacing and normalized to the mean signal intensity of the piriformis muscle to minimize heterogeneity between data sets caused by technical variations in image acquisition.26,28 The IBSI-conform29 and validated software MITK Phenotyping27 were used for radiomics feature calculation. A total of 260 radiomics features were calculated, with 91 first-order features and 169 texture features. Volume and shape features were omitted because these are not expected to carry disease-specific information in this setting, as whole BM spaces are analyzed. Four different radiomics models were trained for the PCI prediction task: once with versus without inclusion of a data set from an older scanner, and once with versus without addition of clinical features. Random forest regression-based models were trained to predict PCI from the radiomics (and clinical) features using the sklearn python package.30 Random forest classifiers based on radiomics (and clinical) features were trained for the predictions of binary variables. All machine learning modeling was performed with Python (Version 3.8.10, Python Software Foundation, Delaware). The models were initially trained on the training set of center 1 only (168 cases) and tested on the center 1 test set (59 cases). Then, the models were retrained on all data from center 1 (227 cases), and finally tested on the external test sets (including a total of 143 cases). Further details on the radiomics analysis are reported in Supplemental Digital Content 4, https://links.lww.com/RLI/A809.

Histological, Cytological, and Cytogenetic Data

Bone marrow biopsies were performed at the posterior iliac crest without image guidance. When assessing PCI from BM biopsy, in line with the International Myeloma Working Group recommendations,14 the higher value from the histological and cytological PCI was used. The 5 cytogenetic aberrations gain(1q), del(13q), del(17p), t(4;14), and t(14;16), which are currently used to form the respective risk stratifications in smoldering MM12 and in MM,10,31 were investigated. In addition to the prediction of the presence of each cytogenetic aberration individually, it was also predicted whether a high-risk cytogenetic status was present. Three different definitions for cytogenetic high-risk were used: high-risk cytogenetic status according to definition 1 (abbreviated “HR-C def 1”), based on the definition of cytogenetic high-risk aberrations in R-ISS [presence of any aberration of the following: del(17p), t(4;14), and t(14;16)]10; “HR-C def 2,” based on the definition of cytogenetic high-risk aberrations in R2-ISS [presence of any aberration of the following: gain(1q), del(17p), t(4;14), and t(14;16)]31; and “HR-C def 3,” based on the definition proposed for smoldering MM [presence of any aberration of the following: gain(1q), del(13q), del(17p), t(4;14), and t(14;16)].12

Statistical Analysis

Dice scores were calculated to quantify the agreement between automatic and manual segmentations. Pearson correlation was used to evaluate the correlation between predicted and actual PCI values. The area under the receiver operating characteristic (AUROC) was used to evaluate the cytogenetic risk status or cytogenetic aberration prediction. Spearman correlation coefficient was used to evaluate the correlation between individual radiomics features and PCI. The Wilcoxon test was used to assess the difference in predicted PCI values between the cytogenetic high-risk and standard risk group, or between the groups with and without the respective cytogenetic aberration. The Gini feature importance was used to report the relative influence of a feature for the prediction model and was calculated as implemented in scikit-learn.30 The 95% confidence intervals (CIs) for Pearson correlation coefficients and AUROCs were calculated. P values <0.05 were considered statistically significant. The statistical analysis was performed with Python (Version 3.8.10; Python Software Foundation, Delaware) and R (Version 4.0.1 R Foundation for Statistical Computing, Vienna, Austria).

RESULTS Study Cohort and Data Sets

A total of 672 MRIs from 512 patients (median, age 61 years; interquartile range, 53–67 years; 307 men) from 8 centers and 370 corresponding BM biopsies were included in this study. An overview over the data sets is displayed in Figure 2. Details on inclusion and exclusion at each stage are reported in the methods and in the supplements (Supplemental Digital Content 1, https://links.lww.com/RLI/A806). Table 1 reports descriptive information for each data set.

TABLE 1 - Description of Study Cohorts A. Data for Segmentation Experiments Training Set for nnU-Net Internal Test Set/Interrater Variability Center 2 Test Set Multicenter Test Set  n MRIs (n patients) 470 MRIs (from 310 patients) 8 wb-MRIs (from 8 patients) 15 wb-MRIs (from 15 patients) 14 wb-MRIs (from 14 patients) Patient characteristic  Male sex, n (%) 186 (60%) 4 (50%) 8 (53%) 8 (57%)  Age in yearsb 61 (53–67) 59 (48–67) 51 (47–66) 60 (54–64) Disease stage  MGUS 26 2 0 0  SMM 131 2 4 0  NDMM 182 4 11 14  ISS I/II/III (n.a.) 114/25/22 (21) 0/2/1 (1) 4/3/4 5/4/5  On/after therapy or n.a. 131 0 0 0 Tumor load surrogates  PCI in %b 23 (12–50; 257) 27 (20–63) 43 (19–59) 48 (15–71)  M-protein in g/Lb 20 (11–35; 130) 30 (15–50) 26 (17–35) 36 (28–44) B. Data for Radiomics Experiments Internal Training Seta Radiomics Internal Test Seta Center 2, High-Quality Test Set Center 2, Other Test Set Multicenter Test Set  n MRIs (n patients) 168 wb-MRIs (from 166 patients) 59 wb-MRIs (from 59 patients) 32 wb-MRIs (from 32 patients) 75 wb-MRIs (from 75 patients) 36 wb-MRIs (from 36 patients) Patient characteristic  Male sex, n (%) 101 (61%) 30 (51%) 16 (50%) 56 (75%) 19 (53%)  Age in yearsb 61 (52–67) 60 (51–69) 64 (53–70) 61 (53–68) 60 (51–65) Disease stage  MGUS 6 4 0 2 0  SMM 47 17 8 11 0  NDMM 115 38 24 62 36  ISS I/II/III (n.a.) 76/18/14 (7) 23/9/2 (4) 13/3/8 (0) 33/13/12 (4) 11/15/10 (0)  On/after therapy or n.a. 0 0 0 0 0 Tumor load surrogates  PCI in %b 23 (12–46) 20 (15–55) 43 (18–68) 26 (10–48) 51 (30–80)  M-protein in g/Lb 24 (12–41; 28) 24 (14–43; 13) 24 (10–36; 7) 32 (15–40; 21) 37 (21–46; 16) Cytogenetics  HR-C def 1c 37% (22/60) 27% (4/15) 23% (7/31) 23% (15/65) 17% (6/35)  HR-C def 2c 56% (40/71) 53% (8/15) 45% (14/31) 47% (31/66) 43% (15/35)  HR-C def 3c 78% (62/80) 74% (14/19) 64% (18/28) 69% (45/65) 49% (17/35)  gain(1q)d 35% (35/101) 33% (7/21) 37% (10/27) 40% (26/65) 35% (12/34)  del(13q)d 48% (49/102) 57% (13/23) 54% (15/28) 58% (38/65) 35% (12/34)  del(17p)d 10% (10/102) 4% (1/24) 7% (2/27) 12% (8/65) 9% (3/35)  t(4;14)d 24% (14/59) 27% (4/15) 19% (5/27) 14% (9/65) 9% (3/35)  t(14;16)d 2% (1/54) 0% (0/16) 4% (1/27) 2% (1/65) 3% (1/35) MRI biopsy interval (days)b 5 (0–29) 17 (2–61) 6 (1–24) 4 (1–13) 5 (2–15)

Table A reports descriptive information for each data set included the segmentation experiments, and Table B reports descriptive information for each data set included in the radiomics experiments.

aNote: Initially, algorithms were trained on only the training set from center 1 and tested on the internal test set. Then, radiomics algorithms were retrained on all data from center 1 and then the resulting models were tested on the external test sets.

bMedian (interquartile range; n missing [only reported in case there were any cases with missing information]).

cPercentage with high risk (number high risk/number all).

dPercentage with cytogenetic aberration (number cytogenetic aberration present/number cytogenetic aberration tested). HR-C def 1, high-risk cytogenetic status according to definition 1 based on the definition of cytogenetic high-risk aberrations in R-ISS [presence of any aberration of the following: del(17p), t(4;14), and t(14;16)]10; “HR-C def 2,” based on the definition of cytogenetic high-risk aberrations in R2-ISS [presence of any aberration of the following: gain(1q), del(17p), t(4;14), and t(14;16)]31; and “HR-C def 3,” based on the definition proposed for smoldering MM [presence of any aberration of the following: gain(1q), del(13q), del(17p), t(4;14), and t(14;16)].12

n.a., not available; MGUS, monoclonal gammopathy of unknown significance; SMM, smoldering multiple myeloma; %, percentage of this cohort; PCI, plasma cell infiltration in the bone marrow in %.

Quality of Automated Pelvic Bone Marrow Segmentation

The Dice scores for the automated pelvic BM segmentations and for the interrater variability between 2 radiologists are reported in Table 2. Figure 3 displays automated segmentations in 5 examples, including cases with severest pathologies.

TABLE 2 - Quality of Automatic Segmentation and Interrater Variability Test Set n nnU-Net vs Radiologist: Right Pelvis nnU-Net vs Radiologist: Left Pelvis Interrater Variabilitya: Right Pelvis Interrater Variabilitya: Left Pelvis Internal test setb 8 0.96 ± 0.03 0.96 ± 0.03 0.88 ± 0.02 0.87 ± 0.02 Center 2 test setc 15 0.93 ± 0.01 0.92 ± 0.01 Multicenter test setd 14 0.89 ± 0.03 0.90 ± 0.03

Mean Dice scores (± standard deviation) are reported to assess quality of automatic bone marrow segmentation, and to compare it with the interrater variability of segmentations between 2 radiologists.

aBetween manual segmentations from 2 different radiologists.

bLast 2 by date for each data set I–IV.

cNewest 15 by date.

dNewest 3 by date per center. From center 4 and center 7, only one processable data set was available. This data set comprises 14 MRIs from 6 centers, acquired with 5 different scanner models from 3 different vendors.

FIGURE 3:

Exemplary automated segmentations in different pathologies and different external centers. On the left, one slice of the MRI scan for each case is displayed. The automatic segmentation (bone marrow of the right pelvis: red, bone marrow of the left pelvis: yellow, medial part of piriformis muscle: blue) is superimposed on the respective image in the middle, and the 3-dimensional model of the automatic segmentations is shown on the right. A, A 68-year-old male patient with normal-appearing, T1w-hyperintense, “fatty” bone marrow. B, A 55-year-old male patient with over 30 T1w-hypointense focal lesions in the right and left pelvis, ranging from 7 mm to 6.2 cm. All focal lesions were correctly included in the segmentation. C, A 66-year-old female patient with severely T1w-hypointense bone marrow. The signal intensity of the bone marrow approximates the signal intensity of muscle, representing a severe diffuse infiltration pattern. This patient showed 70% plasma cell infiltration in the bone marrow biopsy. Despite the severe pathology, the bone marrow was segmented correctly. A–C, Images were all acquired at center 2. D, A 67-year-old female patient with large, secondary extramedullary lesion at the right posterior iliac crest. The scan was acquired at center 5 with a scanner from a different vendor, and consequently, the image appears somewhat different compared with examples A–C. Despite the large focal lesion (6.1 cm by 4.2 cm) with a paramedullary component and the fact that the image was acquired at a scanner from a different vendor, the segmentation was precise and even contained parts of the lesions, which extend beyond the regular shape of the bone. E, A 67-year-old male patient with a large (6.5cm × 4.8 cm), secondary extramedullary focal lesion at the right posterior iliac crest and several focal lesions in the left iliac crest. This scan was acquired at center 6 with an older MRI scanner model, and the image appears somewhat different compared with examples A–C. Despite the considerable paramedullary part of the focal lesion at the right side and the different image characteristic, the automated segmentation was precise and included even the paramedullary parts.

Quantitative Profiling of Bone Marrow Phenotypes Using Radiomics

A wide variety of different morphologic BM patterns can be observed in MRI in patients with monoclonal plasma cell disorders. Figure 4 displays several exemplary cases of these varying MRI BM patterns and the resulting radiomics signature from the pelvic BM for each case. These cases exemplify how differences in morphologic MRI patterns lead to differences in the extracted quantitative, objective radiomics profiles.

FIGURE 4:

Visual bone marrow patterns and resulting radiomics profiles. Panel A displays different visual bone marrow patterns as observed in T1w images, and panel B shows the corresponding radiomics signatures extracted from the pelvic bone marrow for each case. Patient 1 (P1) and P2 show hyperintense, relatively homogeneous bone marrow, which is the physiological pattern observed in elder patients. P3 and P4 show a pattern with T1w-hypointense focal lesions. P3 shows a patient with over 30 focal lesions in the right and left pelvis, sizes ranging from 7 mm to 6.2 cm, representing a pattern of multiple focal lesions with varying size, combined with a homogeneous, intermediate diffuse infiltration. P4 shows a patient with multiple (20) focal lesions in the pelvis, ranging from 5 mm to 1.6 cm, and a patchy diffuse infiltration. P5 and P6 show 2 cases with severely T1w hypointense bone marrow, which approximates the signal intensity of muscle tissue, representing a severe diffuse infiltration. In line with this imaging finding, patients showed 70% (P5) and 85% (P6) plasma cell infiltration in the bone marrow. P7 to P10 represent subtypes of intermediate diffuse infiltration, which would be neither classified as normal, focal, nor severe diffuse infiltration. P7 shows a pattern with very small nodular components in which most nodules have relatively similar size, and a relatively low contrast between the hypointense and hyperintense nodules. P8, in contrast, shows a pattern in which fewer nodules are present, and nodules are on average larger and show greater variation regarding size. In addition, in contrast to P7, at least some of the nodules are more hypointense (approximating signal intensity of muscle), and the pattern is characterized by a higher contrast between the hypointense nodules and the fatty bone marrow. P9 shows a case with a relatively homogeneous, intermediate T1w signal intensity, which lies in between the signal intensity in normal fatty bone marrow and severe diffuse infiltration. P10 shows a case in which homogeneous fatty plates and interspersed, confluent hypointense insulae are observed, in which the insulae are not as hypointense as muscle.

The differences between MRI morphologically normal-appearing BM, focal lesion pattern, and severe diffuse infiltration are unequivocal (Fig. 4, P1–P6). Besides those, there is a large group of patterns that might be classified as intermediate diffuse infiltration (Fig. 4, P7–10). These are clearly heterogeneous and therefore might represent different MM disease subtypes. However, the complexity of such patterns can hardly be reported in a structured, reproducible manner based on a visual assessment by radiologists. Detailed, objective analysis of these complex patterns and their systematic correlation with tumor tissue characteristics (assessed by BM biopsy at the iliac crest) are an optimal use case for machine learning–based image analysis approaches as deep learning and radiomics.

Automatic Prediction of Plasma Cell Infiltration

Four different models were trained and tested. Once, only data sets I to III of the internal data sets were used, whereas data set IV, which was acquired with an older scanner and had markedly lower image quality, was omitted. Second, all data sets I to IV from the internal data set were included. In both scenarios, one model was trained on radiomics features only, and an additional model was trained both on the radiomics features and the clinical features age and body mass index. The correlations between the predicted PCI values and the actual PCI values for the different prediction models are reported in Table 3. The model based on data set I–III using only radiomics features without clinical parameters showed the highest correlation coefficient between predicted PCI and actual PCI (r = 0.71, P < 0.001) on the internal data set. The correlation between predicted PCI and actual PCI on the external data sets was worse than in the internal data set. However, on all external data sets, the model trained on radiomics features from data set I–III without clinical parameters predicted PCI values which were significantly correlated to the actual PCI values (all P's ≤ 0.01), with correlation coefficients between 0.30 and 0.56. The model including additional clinical features performed quite similar to the model without clinical features, with the main difference that in the data set from center 2 with variable imaging quality, it performed somewhat better (r = 0.38 vs r = 0.30). Addition of data set IV to enhance the training data set did not markedly change the performance of the PCI prediction models, neither when using only radiomics features nor when additionally including the clinical features age and body mass index. As a benchmark for the interpretation of the correlation coefficients between predicted PCI and actual PCI, we investigated the correlation between PCI from histological and cytological assessment, and found the correlation coefficient to be 0.53 (P < 0.001).

TABLE 3 - Accuracy of the Prediction of Plasma Cell Infiltration of Different Models Test Set Model 1: Trained on Data Set I–III; Radiomics Features Model 2: Trained on Data Set I–III; Radiomics and Clinical Features* Model 3: Trained on Data Set I–IV; Radiomics Features Model 4: Trained on Data Set I–IV; Radiomics and Clinical Features* Internal   r 0.71 [0.51, 0.83] 0.66 [0.44, 0.80] 0.56 [0.35, 0.71] 0.56 [0.35, 0.72]   p <0.001 <0.001 <0.001 <0.001 Center 2, high-quality subset   r 0.45 [0.12, 0.69] 0.42 [0.08, 0.67] 0.42 [0.09, 0.67] 0.38 [0.04, 0.64]   p 0.009 0.02 0.02 0.03 Center 2, other subset   r 0.30 [0.07, 0.49] 0.38 [0.13, 0.59] 0.22 [−0.01, 0.43] 0.39 [0.15, 0.60]   p 0.01 0.004 0.06 0.003 Multicenter, test set   r 0.57 [0.30, 0.76] 0.58 [0.07, 0.85] 0.45 [0.15, 0.68] 0.58 [0.08, 0.85]   p <0.001 0.03 0.006 0.03

The correlation coefficient r and the P values from Pearson correlation are reported for each prediction model on each test set.

*Composition of test sets center 2, other and multicenter deviated from the other analyses, as body mass index was not available for all patients in these data sets.

Figure 5 visualizes the 15 most important radiomics features (according to the PCI prediction model based on the internal training data set using radiomics features only) for the internal training set and the internal test set and provides a quantitative analysis how each of the features is correlated to the PCI in each data set. When investigating visual patterns in the radiomics heat map of the internal training set, no general, continuous trend of each radiomics feature from patients with low PCI toward patients with high PCI can be observed. Rather, especially in patients with low to intermediate PCI, different patterns in the radiomics signatures can be found. The fact that quite heterogeneous quantitative imaging features are found in patients with low to intermediate PCI is very much in line with the observation that there are very heterogeneous visual patterns in the BM in such patients, as demonstrated in the exemplary cases in Figure 4. However, especially in patients with (very) high PCI, a rather distinct pattern in the radiomics heat map becomes apparent (Fig. 5). This comprises lower features values for features as “first-order numeric mode value”/“first-order histogram mode value,” “first-order numeric 30th percentile,” or “first-order numeric minimum,” representing a rather hypointense BM signal.