A retrospective study of MRI scans performed at diagnosis and follow-up in the Unit of Lysosomal Disorders and evaluated by the Spanish Group of GD from April 1995 to May 2022 was conducted. A total of 441 bone marrow MRI examinations (S-MRI protocol) were included from 131 patients diagnosed with GD. Infiltration in the lumbar spine, pelvis, and femora was evaluated, according to the signal intensity on T1 and T2 WI. Progressive values for each MRI pattern were assigned following S-MRI description: homogeneous (H, 4 points), non-homogeneous diffuse (NHD, 3 points), non-homogeneous mottled (NHM, 2 points), or non-homogeneous reticular (NHR,1 point), and normal or no infiltration (N,0) (12). The existence of complications (infarcts, necrosis, fractures, arthropathy, or bone crisis) was also taken into account (4 points). We applied the structured bone marrow MRI report in each imaging study (16). The report model and S-MRI description are included in Supplementary Figure 1.

All studies were reevaluated by the same expert radiologist in a blinded fashion to ensure that the structured report model and analysis were conducted as objectively as possible.

We also analyzed the information reported in each patient's clinical record, collecting demographic, genetic, clinical, and analytical data; the type of treatment; and the accumulated years of treatment exposure. The variables analyzed are listed in Supplementary Table 1.

To assess the evolution of bone involvement over time, the studies were divided into four groups. The first group (group A) included studies performed at diagnosis or before starting therapy (baseline). The second group included studies performed between the 1st and 4th year of follow-up (group B). The third group included studies performed between the 5th and 9th year of follow-up (group C), and the fourth group included studies performed after 10 years or more of follow-up (group D).

As part of the evaluation of bone involvement, bone mineral density was also estimated considering the criteria recommended by the WHO for the diagnosis of osteopenia and osteoporosis. The T-score or Z-score was used accordingly (17). Gender, genotype, spleen status, type of therapy, age at the start of treatment, and accumulated years of treatment in different subgroups were also considered.

A descriptive study was carried out, presenting qualitative variables as percentages. Mean and standard deviation were presented for those quantitative variables that followed a normal distribution, and median, interquartile range, or range between maximums and minimums were presented for those that did not follow a normal distribution. In addition, a correlation analysis was undertaken between numerical variables using Pearson's linear correlation index, and for categorical variables, the X² test was performed.

The data were analyzed using the IBM SPSS statistics 27.0 statistical program.

For the comparative analysis of the results, the Student's t-test for the comparison of independent samples was used for quantitative variables. To determine the suitability of this test, an analysis of the normality of the distributions was first performed.

From multiple random forest test models, three were selected and trained to identify features that can predict the risk of bone complications. Bone complications were defined by the presence of intraosseous ischemic events (bone crisis, infarcts, avascular necrosis, and fractures) during the follow-up. Model A included all variables that were described as significant in a previously published study from our group using the same population (18), model B considered whether a treatment was applied or not, and model C ignored the S-MRI punctuation score.

The model parameters were optimized using a grid search. As random forest modeling uses bootstrapping, no cross-validation was used to reduce the training time. Model performance on the validation dataset was evaluated using the ROC curve, AUC, accuracy, and f1-score. The model was created using the scikit-learn package for Python 3.10.4 (19).

The cohort included 131 GD patients (62 female patients and 69 male patients) and 310 follow-up bone marrow MRI studies. The general characteristics of the patients are detailed in Table 1.

At the time of the first MRI, the median S-MRI was 8.4 (95% CI 0–25) (Figure 1C). The median S-MRI at baseline according to age groups was years: 0–20: 8.78; 21–40: 7.78; 41–60: 7.68; >61: 11.0; S-MRI according to gender distribution and age group is shown in Figure 1A.

Regarding gender distribution, the median S-MRI at baseline was higher in male patients [9.10 (95% CI 0–25)] than in female patients [7.71 (95% CI 0–24)], and this difference was significant (p < 0.001). During the follow-up, the reduction in bone marrow infiltration was different between men and women. A 20.8 % decrease in S-MRI was observed in male patients after 10 years on therapy, while female patients achieved a 39.0% reduction in S-MRI earlier (5–9 years) (Figure 1B). No significant differences in genotype distribution, age, and spleen status were observed between genders (see Table 1). Globally, the maximum reduction was observed after 5–9 years on therapy, and it remained stable after 10 years of follow-up (Figure 1C).

Bone marrow infiltration in all locations was found in 80 (61.5%) patients: 38 (47.5%) were female patients, and 42 (52.5%) were male patients. The MRI pattern was homogeneous in the lumbar spine in 47 (36.1%) patients, in the pelvis in 14 (10.7%) patients, and in the femurs in 5 (3.8%) patients (Figure 2A). At baseline, 63 patients (48.4%) showed complications such as bone crisis, fractures, AVN, or infarcts.

According to genotypes, patients were classified into homozygous c.1226A>G (N370S/N370S), heterozygous c.1226A>G/c.1448T>C (N370S/L444P), heterozygous c.1226A>G/other variant, and other different variants (Table 1). The analysis at baseline of the degree of BM infiltration and bone complications according to genotype showed no significant differences between 13 homozygous c.1226A>G (N370S/N370S) and 44 heterozygous c.1226A>G/c.1448T>C (N370S/L444P), but a significant difference between the 57 heterozygous c.1226A>G/other and 17 other/other genotype patients with different variants (p = 0.017) was observed (see Supplementary Figure 2).

There were 310 follow-up studies of these patients available for evaluation, which were analyzed according to the specified time groups. The patterns reported at baseline in the different locations (S: spine; P: pelvis; and F: femur) are represented with color code in Figure 2A. The changes (%) in the infiltration pattern during follow-up at each site using the same color code according to the groups were A (baseline), B (1–4 years follow-up), C (5–9 years follow-up), and D (10 or more years follow-up) (Figure 2B).

The analysis of changes in infiltration patterns related to the type of treatment was difficult due to the diversity of treatments used over time and the diversity of applied doses of enzyme replacement therapy. To obtain some useful information, we classified the patients generically in enzyme replacement therapy (ERT) distributed in their follow-up groups and we differentiated between patients treated with miglustat and those treated with eliglustat (Table 1).

The median S-MRI in 31 patients with spleen removal at baseline was 13.16 (95% CI 0–25), which was significantly higher than that in non-splenectomized patients (6.96, 95% CI 0–24) (p < 0.001). At follow-up, the differences persisted with a median S-MRI in group D (10 or more years) in non-splenectomized patients of 5.42 (95% CI 3.8–6.3) vs. 9.02 (95% CI 4.5–12.7) in splenectomized patients (p < 0.001). Patients with spleen removal also had a significantly higher incidence of bone complications (40.8 vs. 16.0%; p 0.0001). No patients were splenectomized during follow-up.

Bone mineral density was evaluated in 119 (90.8%) patients following the WHO criteria using the Tor Z-score when applicable (15). In the first study, osteoporosis was defined in 28 (23.5%), osteopenia in 25 (21.0%), and was normal in the rest (66; 55.5%). A decrease in BMD was observed in 74.6% of splenectomized patients compared to 43.8% of non-splenectomized patients (p < 0.001).

Bone mineral density data were available from 256 follow-up studies, of which 170 corresponded to patients under 50 years of age and 86 to those over 51 years old. In terms of gender, 45.0% of male patients and 56.4% of female patients under 50 years presented a decrease in BMD, and this difference was not significant (p = 0.4). In the group over 51 years of age, there was a decrease in bone mineral density of 44.0% in male patients and 61.6% in female patients; this difference was significant (p = 0.04). Analysis of the follow-up of BMD in the different time periods showed no changes in the male patients, while in female patients there was a significant increase in the loss of BMD throughout the follow-up period in this study (p = 0.001).

There was no evidence of a relationship between decreased BMD and the presence of intraosseous vascular events during follow-up.

Analysis of the GBA1 genotype within the different subgroups showed that patients homozygous for c.1226A>G (N370S) had a significantly lower prevalence of lost BMD compared to the rest of the genetic subgroups (p = 0.003).

The random forest model A achieved an AUC of 75.82% with an accuracy of 78.10% and an f1-score of 75.18%. We obtained a true positive rate (TPR) of 69.70% and a false positive rate (FPR) of 18.06% with a decision boundary threshold of 0.5 (Figure 3A). Using the mean decrease in accuracy as the feature importance metric, the most important features for this model were the S-MRI, the age at first treatment, and the treatment used. The random forest model B, which considered whether a treatment was applied or not, achieved an AUC of 85.73% with an accuracy of 83.81% and an f1-score of 87.21%. We obtained a true positive rate (TPR) of 90.91% and a false positive rate (FPR) of 19.44% with a decision boundary threshold of 0.5 (Figure 3B). Using the mean decrease in accuracy as the feature importance metric, the most important features for this model were the S-MRI, the age at first treatment, and the extent of spine infiltration.

In order to understand the importance of each variable in the developed model, new models were generated removing one single variable at a time. The variable with more relevance in the accuracy was the S-MRI. Model C did not contain the S-MRI score and had a substantial drop in accuracy; it achieved an AUC of 69.76% with an accuracy of 74.29% and an f1-score of 69.92%. We obtained a true positive rate (TPR) of 57.58% and a false positive rate (FPR) of 18.06% with a decision boundary threshold of 0.5 (Figure 3C). Using the mean decrease in accuracy as the feature importance metric, the most important features for this model were the femur infiltration, age at first treatment, and the diagnostic age.

For almost 30 years, GD patients have been treated with ERT, which is effective in clearing substrate deposits, especially at the visceral level. However, bone involvement and its complications (bone infarcts, avascular necrosis, fractures, bone crisis, cortical thinning, osteopenia, and osteoporosis or lytic bone lesions) are problems that have not been completely solved with current treatments (2, 20, 21). The progressive storage of glucocerebroside in the bone marrow, cytokine imbalances, and vascular compromise are some of the proposed explanations for the development of bone complications, although there is no clear explanation about why it appears in some patients and does not in others with similar clinical and genetic characteristics (22). Bone complications are the major cause of morbidity and one of the most debilitating aspects of GD (1, 3, 21–23). Despite various efforts to study different biomarkers to predict the development and intensity of bone involvement and its complications, no specific markers have been identified (24–28).

At baseline, our results are in accordance with previous publications that show severe bone involvement in splenectomized patients (18, 29). In addition, these patients had a significant loss of BMD, which may increase the risk of fractures in this subgroup of patients although this estimate should be taken with caution as there may be a bias in the estimation of BMD due to the presence of infarcts, especially in the femoral neck (30). Interestingly, there was a significant difference in the S-MRI at baseline between male (9.10) and female patients (7.70) (p = 0.047), which was not previously reported. This could be justified in part by the differences in bone maturation according to sex (31), which have been reported in MRI imaging performed in healthy subjects, mainly in the lumbar spine. MRI imaging showed that male and female subjects convert hematopoietic marrow to fatty marrow in the lumbar vertebral bodies in significantly different ways (32). Nevertheless, other factors, such as hormonal or vascular changes, may be involved. In addition, clinical characteristics (spleen status, genotype, and age) may have some significance; however, in our analysis, non-significance was found. Other studies were carried out with different methodologies such as the Dixon quantitative chemical shift imaging (QCSI reported by M Maas et al. to assess bone marrow in the spine of healthy subjects) (33). This technique allows assessing the fat fraction of the bone marrow with greater accuracy applied to a normal population stratified by age/gender; the results are in agreement with those previously reported by Ishijima et al. (33). The fat fractions in female subjects slowly increases with age until 44 years and rapidly after 45 years. In males, there is a rapid increase in the fat fraction up to the age of 25 years and it stabilizes in the following decades up to the age of 60 years. More recently, another tool has been described to measure the fat fraction in the lumbar spine. Fat fraction quantification of the bone marrow in the lumbar spine using the LiverLab assessment tool with results superimposable to those obtained using QCSI (34).

Our follow-up data show that the maximum reduction in BM infiltration occurs between 5 and 9 years on therapy, subsequent studies showed stability, and this is in line with other groups' observations (35). As expected, low baseline infiltration in women was associated with faster clearance.

No significant differences were identified regarding the age at diagnosis, age at first therapy, and the number of complications between male and female patients. However, the differences according to gender persisted during follow-up, with a median S-MRI after 10 or more years of 7.29 in male patients vs. 6.60 in female patients (p < 0.001), this was independent of bone mineral density status (1, 10, 18).

There is enormous variability in bone marrow patterns between age groups; some conditions, such as post-bleeding anemia or therapies, may complicate bone marrow image interpretation (36–38). In any case, training is required and uniformity in the description of its assessment by the radiologist is desirable (12). The application of the structured report template in our cohort for BM MRI studies improved standardization and the quality of the radiology reports, allowing easy comparison and the incorporation of data into the machine-learning system.

In the area of rare diseases, the use of machine learning provides an opportunity to analyze agglomerated and heterogeneous data to create quality predictive models and identify risk features (39). This can be useful to improve the study of small cohorts of patients (40–42) and facilitate differential diagnosis; recently, its application in neuromuscular diseases was reported (43).

In line with our previous study (18), the best-generated model (model A) identified the S-MRI score and the age at first treatment as risk predictors for developing advanced bone disease. In this model, femur and pelvis infiltration at baseline and homogeneous pattern infiltration at any location also predicted a greater degree of bone involvement. It is logical that infiltration in locations, such as the pelvis and femur, is related to the extent of the disease. In addition, as expected, the intensity of the infiltration pattern (homogeneous > non-homogeneous) had an impact on the severity of bone involvement. These aspects are not considered independently at the initial calculation of the overall staging of the disease (5) but they are included in the S-MRI score (10, 12). However, they could be potential independent risk features for developing advanced bone disease and support the decision of early treatment (44). Further validation in another cohort is warranted.

We have observed that S-MRI is consistently the most important variable across all the models developed. Age and treatment contributed to the final accuracy indeed but to a lower extent, and their contribution is not consistent in all of the developed models. Model A was demonstrably the best of the three and includes information about all the different therapies. Model A could be therefore implemented into a webpage algorithm or a service for clinicians to feed it with the data and obtain a potential prediction of the prognosis of the patient.

In addition, this study supports the importance of early diagnosis to allow for tailored therapy. Delays in the diagnosis and therapy were related to bone complications in this cohort.