A validated disease severity scoring system predicts risk of fractures in treated adults with type 1 Gaucher disease

INTRODUCTION

Gaucher disease type 1 (GD1; OMIM #230800) is an autosomal recessive lysosomal storage disorder caused by deficient activity of glucocerebrosidase due to pathogenic variants in the GBA gene (GBA1; G/L MIM 606463^[1]. GD1 exhibits marked clinical heterogeneity, ranging from asymptomatic disease to progressive involvement of hematologic, visceral, and skeletal systems. Enzyme replacement therapy (ERT) has transformed the natural history of GD1; however, responses - particularly with respect to skeletal outcomes - remain variable^[2-4].

The skeletal manifestations of GD1, chronic bone pain, acute bone crisis, osteonecrosis and fragility fractures contribute substantially to morbidity and reduced quality of life^[5,6]. Fragility fractures are related to skeletal pathophysiology including abnormalities of the bone marrow and the effects of infiltration with glucocerebroside-engorged macrophages (Gaucher cells); abnormalities of trabecular and cortical bone including modeling and remodeling defects; loss of bone mineral mass with osteopenia/osteoporosis or lytic lesions; vasculopathy and osteonecrosis that adversely affects bone structural and mechanical integrity^[2,7,8]. Displacement or suppression of bone marrow adipocytes may also inhibit osteoblastic synthesis of new bone shifting the equilibrium in favor of bone resorption and bone mineral loss^[2].

Estimation of fracture prevalence is challenging due to disease and event rarity, heterogeneous reporting, and limitations of GD1 registry data. Validity of comparisons with control groups is challenged by variable geographical norms, and changing demographics including age distributions, racial and ethnic composition, and imbalances in economic and nutritional conditions^[9]. Retrospective registry studies are compromised by case under-reporting, cross-sectional rather than longitudinal reporting, incomplete data, and inclusion of fragility and non-fragility fractures^[10]. Notwithstanding these caveats, among 706 never treated or pre-ERT patients with GD1 enrolled in the International Collaborative Gaucher Group (ICGG) Gaucher Registry before 1998, 15% had a history of fracture^[11]. In a multi-center study in the United Kingdom (UK), 28% reported fragility fractures^[12]. Among 105 adult GD1 patients in France, 18% had non-vertebral fractures and 15% had vertebral fractures, and in a more recent study of 239 patients, there were more pathological fractures but fewer vertebral compression fractures {58 [24.3%], 23 [9.6%]} within the first 15 years after starting ERT^[13,14]. In comparison, from 2017-2020, in US individuals 50 years or older, the prevalence of fragility fractures (hip, wrist, vertebral) was only 0.8%, 2.2%, 0.4% and non-fragility fractures 0.6%, 1.8%, 0.9%. For individuals younger than 50 years of age, the prevalences for all fractures were 0.7%, 9.8%, 3.0%^[15]. This disparity confirms that at all ages, individuals with GD1 have an inordinate number of fractures compared to those non-affected.

Fracture risk assessment in GD1 has relied largely on dual-energy X-ray absorptiometry (DXA), despite the complex skeletal pathology characteristic of this disease^[16]. In a case-controlled study in which 319 GD1 patients with a first-time fracture were matched for gender, year of birth, treatment status, and splenectomy status with 1,233 GD1 patients without a fracture history, only lumbar spine Z-scores ≤ -1 were identified as a risk factor for fracture [Objective Response (OR) 5.55]^[17]. Nevertheless, lumbar spine Z score as a standalone parameter for assessing fracture risk in GD1 is problematic. In children with GD1, DXA bone density measurements may be underestimated due to low bone size and growth retardation. In adults older than 50 years, T scores rather than Z scores are a more accurate indicator of osteopenia/osteoporosis and fracture risk. DXA scores may be falsely elevated in the presence of osteophyte formation, bone sclerosis, disk space narrowing, and spondylolisthesis. Consequently, the Gaucher Risk Assessment for Fracture (GRAF) score was developed to estimate fracture risk in patients starting ERT or Substrate Reduction Therapy (SRT) based on demographic and disease-related factors including sex, age at treatment initiation, time interval between age at diagnosis and age at treatment initiation, splenectomy status, and pre-treatment history of bone crises [Figure 1]^[18]. GRAF score is used most effectively only in conjunction with DXA and other non-GD-specific risk factors (e.g., smoking, alcohol intake, frailty).

Figure 1. Weighted calculation of DS3 scores for the bone, hematologic, and visceral domains. Note the change in the assessment of bone mineral density as explained in the text. Reproduced from Weinreb et al.^[19]. DS3: Disease severity scoring system; MRI: magnetic resonance imaging; DXA: dual-energy X-ray absorptiometry; BMB: bone marrow burden; CT: computed tomography; AVN: avascular necrosis.

Separately, the Disease Severity Scoring System (DS3) integrates multisystem manifestations of GD1, including skeletal involvement and DXA, and has been validated as a measure of overall disease burden and treatment response^[18,19]. Whether DS3 predicts fracture risk, and how it compares with or complements GRAF has not been previously examined. Therefore, using retrospective data from the DS3 validation study^[19], we investigated whether pre-treatment DS3 scores predict fracture risk in adults with GD1 receiving ERT, whether DS3 outperforms GRAF, and whether combining both measures improves predictive accuracy.

METHODS

This study was performed as part of the research titled “Validation Study of a Severity Score Index (DS3) for Adult Patients with type 1 Gaucher disease” (URFLSD-.2010-01) after review of the protocol (WIRB® Protocol #20100371) and approval of the patient informed consent form dated March 23, 2010 by the Western Institutional Review Board (WIRB®, now WCG IRB Connexus, Seattle, WA, USA). All consent forms were signed by the patients. The signed consent form includes patient permission to publish aggregate anonymized results in scientific journals.

Patient recruitment, eligibility, and calculation of the DS3 scores have been described previously^[19]. The DXA and Magnetic Resonance Imaging (MRI) bone marrow infiltration lack of availability was anticipated when the DS3 score was initially formulated and the score, when appropriately modified, was validated even in the absence of these two elements^[18]. Patient eligibility was contingent on enrollment in the ICGG Registry between 1991-2009. IRB-approved consent forms allowed future retrospective investigations. The study cohort comprised 133 adult patients with GD1. Demographic characteristics are shown in Table 1. Comparison with other large population ICGG studies^[16,20] indicates that the DS3 study population is demographically representative in terms of sex, age at diagnosis, age at first treatment, GBA1 genotype distribution, and percentage of patients with a history of splenectomy [Table 2]. DS3 scores were calculated retrospectively using registry and study-specific data according to published criteria [Figure 1 and Supplementary Table 1]^[18]. Disease severity was categorized as mild (< 3), moderate (3-5.9), or marked (≥ 6). Fractures are entered in the Registry database based on radiological evidence as reported by the site investigator. A distinction among pathological, fragility, or non-fragility fractures is not included because detailed clinical information about the circumstances is not captured in the case report forms.

GRAF scores were calculated for each of the 133 patients enrolled in the DS3 study with the tool illustrated in Figure 2 and per the instructions in Deegan^[16].

Figure 2. Algorithm for calculation of the Gaucher Risk Assessment for Fracture (GRAF) score. Reproduced with permission from Deegan et al.^[16].

Statistical analyses included non-parametric tests, χ² testing with Holm correction for multiple comparisons, Kaplan-Meier survival analysis, and sensitivity analyses. In fracture-free Kaplan-Meier (KM) survival analyses, statistical significance was defined as a two-sided P value < 0.05 after Holm correction.

Receiver Operating Characteristic (ROC) analysis was used to assess the ability of baseline disease severity indices to predict fractures occurring during ERT in our GD1 patient population. Three predictors were evaluated: (1) the continuous DS3 score, (2) the ordinal GRAF score (1 = mild, 2 = moderate, 3 = severe), and (3) a combined model incorporating both measures via multivariate logistic regression.

ROC curves were constructed for each predictor and stratified by age at treatment initiation (< 18, 18-49, and ≥ 50 years). The area under the ROC curve [area under the curve (AUC)] was used as a summary measure of discriminative performance. The optimal threshold for classification was determined using Youden’s index (J = sensitivity + specificity - 1). To assess statistical uncertainty, 95% confidence intervals (CIs) for AUC values were derived using nonparametric bootstrapping (2,000 resamples). For comparison between predictors, a bootstrap-based approximation of the DeLong test was used to estimate differences in AUC (ΔAUC) between the combined model and each single predictor, with one-sided p-values indicating the likelihood that the observed improvement occurred by chance. AUCs ≥ 0.70 were interpreted as moderate discrimination, while AUCs ≥ 0.85 were considered strong discrimination.

Time to first fracture was also analyzed using Cox proportional hazards regression, with follow-up calculated from baseline assessment to first fracture or last follow-up. Age was entered as a continuous variable because of its established association with fracture risk in the general population. Separate multivariable models were constructed to evaluate baseline DS3 and GRAF in relation to fracture risk after adjustment for age. An exploratory model including age, DS3, and GRAF simultaneously was also examined. Hazard ratios (HRs) and 95%CIs were estimated. Additional sensitivity analyses further adjusted for sex, prior fracture history, and splenectomy status.

All analyses were conducted in Python (v3.11).

RESULTS

One hundred thirty-three treated patients were enrolled. 132 were treated with alglucerase or imiglucerase ERT and 1 with miglustat. The miglustat-treated patient was included in the demographic analyses but was excluded from the fracture-free survival analyses which were restricted to patients treated with ERT. When patients were started on treatment, neither velaglucerase alfa, taliglucerase alfa, nor eliglustat were commercially available although some were switched to velaglucerase alfa between 2009 and 2012.

Of the 133 evaluable patients, 33 (24.8%) experienced at least one fracture during follow-up. Patients with fractures were significantly older at ERT initiation, had longer delays between diagnosis and treatment, and had higher baseline age-adjusted Charlson co-morbidity scores^[21,22] and higher DS3 and GRAF scores than those without fractures [Table 3]. Women and patients with prior splenectomy were overrepresented in the fracture cohort. However, there were no inter-group differences regarding GBA1 genotype distributions or baseline weight or body mass index (BMI).

Two-thirds of the fractures involved the spine, hip, or femur [Table 4]. Twelve patients experienced subsequent fractures, typically at older ages (half ≥ 67 years old). Fracture distribution was similar to that reported in larger registry cohorts^[16].

Baseline DXA data were incomplete in a substantial proportion of patients, most of whom enrolled in the ICGG Registry between 1991-2000 when DXA often was unavailable. This prevented any assessment of the value of pre-treatment lumbar spine Bone Mineral Density (BMD) for predicting overall fracture risk in our patient cohort^[17]. Follow-up DXA measurements (osteopenia-osteoporosis combined per Khan et al.^[17]) did not differ significantly between fracture cohorts (~ one year prior to a fracture event) and non-fracture cohorts (last follow-up date) (P = 0.85), indicating that as recently shown by Kahwati et al.^[23], low Aerial Bone Mineral Density (aBMD) alone was insufficient to explain fracture occurrence [Table 5].

Furthermore, after institution of ERT, improvement in DS3 from pre-treatment scores was disassociated from fracture occurrence [Figure 3]. 23 of 33 fracture patients (69.7%) had pre-treatment DS3 scores in the severe range (DS3 ≥ 6). Five of these had fractures during the first year of treatment. Of the 28 patients whose fractures occurred more than one year after starting ERT, 13 (46.4%) had minimal clinically important improvement in DS3 score as defined in Weinreb^[19], before suffering a fracture. Among 10 patients with DS3 < 6, one (baseline DS3 3.25) had a fracture during the first year of ERT (DS3 0.8). Thus, as with DXA, a clinically significant improvement in DS3 did not prevent fracture events in nearly half of the patients who started treatment with DS3 scores indicative of high severity GD1. However, the pre-treatment DS3 score is strongly predictive of continued risk for fractures even after commencement and continued ERT.

Figure 3. Changes in DS3 scores from pre-treatment (baseline) to one year prior to fracture occurrence. The sub-population in this figure comprises the 33 patients from the total study population of 133 patients who had a fracture after initiation of treatment for GD1. The baseline DS3 scores are shown in blue. DS3 scores recorded approximately one year before a fracture occurred are shown in red. Five patients (designated with an “X” below the patient number) had a fracture during the first year of treatment and were excluded from the subsequent analysis of change in DS3 score. Of the 28 remaining patients, 13 (46.4%) had a decrease in DS3 ≥ 3.15, indicating a minimal clinically important improvement as defined in the literature^x, despite which they nonetheless suffered a fracture. These thirteen patients are designated with asterisks. The vertical green line separates the 23 patients with a baseline DS3 ≥ 6.00 (severe severity) from the 10 patients with a lower baseline score. Of these, one patient had DS3 < 3.00. The other nine had baseline DS3 3.00-5.99 (moderate severity). DS3: Disease severity scoring system.

Fracture-free survival differed significantly across DS3 severity category [Figure 4 and Supplementary Figure 1], with higher fracture risk among patients with DS3 ≥ 6. This result supports our hypothesis that DS3 that includes bone variables such as DXA, chronic bone pain and crises, osteonecrosis, lytic lesions, and qualitative severity of bone marrow infiltration with Gaucher cells, might outperform GRAF for predicting future fracture risk in treated patients with GD1. However, inclusion of history of pre-treatment fractures could cause confounding due to predictive circularity. Thus, we performed a Kaplan-Meier sensitivity analysis in which fracture-free survival in all patients with a history of pre-treatment fractures (N = 32) was compared to all those with no history of fracture (N = 100). There was no difference between the groups (log rank P = 0.9; Figure 5A). Furthermore, a similar analysis restricted to 58 patients with a baseline DS3 ≥ 6 in whom the on-ERT fracture risk at 10 years was 30% [Figure 4], also showed no difference between those who had a history of pre-ERT fracture and those who did not [Figure 5B]. Thus, it is unlikely that history of prior fracture invalidates DS3 as a predictor of future fracture risk. Furthermore, in exploratory sensitivity analyses associated with the multivariable Cox model that adjusted for sex, prior fracture, and splenectomy status, DS3 remained associated with fracture risk.

Figure 4. Kaplan-Meier plot of fracture-free survival stratified by DS3 score category. DS3 < 3 (mild); DS3 3- < 6 (moderate); DS3 ≥ 6 (severe). The number of patients at risk at the start of study and after 5, 10, 15, 20 years is shown above the X-axis. P values were calculated by Holm-corrected log rank-χ² tests with results including the hazard ratios, 95% confidence intervals, and the Holm P values shown in the figure. Based on the Holm P values, the fracture-free survival for the tertile with the highest DS3 scores is significantly less than for the moderate and mild DS3 severity tertiles. The difference in fracture-free survival between the moderate and mild DS3 severity tertiles is borderline, but not statistically significant (Holm P = 0.069). DS3: Disease severity scoring system; ERT: enzyme replacement therapy; CI: confidence interval; HR: hazard ratio.

Figure 5. Sensitivity analysis. (A) Kaplan-Meier fracture-free survival plots for 32 patients who had a history of fracture before starting treatment with 100 patients who had no history of pre-treatment fracture. The number of patients at risk at the start of study and after 5, 10, 15, 20 years are shown above the X-axis. P values were calculated by Holm-corrected log rank-χ² tests; (B: post-treatment fracture-free survival in 25 patients with a history of pre-ERT fractures and 32 with no such history. This analysis was restricted to patients with baseline DS3 ≥ 6 (i.e. severe disease). DS3: Disease severity scoring system; ERT: enzyme replacement therapy; HR: hazard ratio.

In the GRAF score Kaplan-Meier fracture-free survival analysis [Figure 6], although the overall log-rank test was significant, pairwise comparisons among GRAF tertiles were not statistically significant after Holm correction. Although the fracture-free survival curves suggest that DS3 outperforms GRAF for fracture risk assessment, the fracture risks at 5, 10 and 15 years are not that discrepant [Table 6]. The differences in the fracture-free survival curves may partly be a function of patient population differences between the DS3 and GRAF sub-groups as each GRAF tertile had an admixture of patients with variable DS3 scores and vice versa [Figure 7]. More importantly, the GRAF score calculation is age dependent and, unlike DS3, patients are assigned a different risk category according to the age at which treatment for GD is initiated (< 18 years; 18-50 years; > 50 years)^[16]. This accords with accepted knowledge that, increasing chronological age, especially after age 50-60 years, is an important fracture risk factor, although it is a complex surrogate biomarker comprising menopausal status in women, prevalence of falling, nutritional status, frailty, level of physical activity, and concurrent illnesses and medications ^[24].

Figure 6. Kaplan-Meier plot of fracture-free survival stratified by GRAF score risk category. 3 = marked; 2 = moderate; 1 = mild. The number of patients at risk at the start of study and after 5, 10, 15, 20 years is shown beneath the X-axis. P values were calculated by log rank-χ² tests as explained in the legend to Figure 4. P values for pairwise comparisons are Holm-corrected; although the overall log-rank test is significant, pairwise differences between groups were not statistically significant after correction. GRAF: Gaucher risk assessment for fracture; ERT: enzyme replacement therapy; CI: confidence interval; HR: hazard ratio.

Figure 7. (A) Distribution of GRAF scores by DS3 category (3,2,1); (B) Distribution of DS3 severity groups by GRAF score risk category (3,2,1). For GRAF risk categories, 3 = high; 2 = intermediate; 1 = low. For DS3 severity categories, 3 = marked; 2 = moderate; 1 = mild. GRAF: Gaucher risk assessment for fracture; DS3: disease severity scoring system; GD: Gaucher disease.

We therefore separated each DS3 severity group (severe, moderate, mild) and assessed the risk of fracture 10 years after start of ERT depending on whether treatment was initiated < 18 years, between 18-49 years, or ≥ 50 years [Table 7]. The only patients at risk for fracture were those older than 50 years with baseline DS3 scores indicative of moderate or high severity and younger adults with baseline DS3 scores ≥ 6. Forty-nine percent of the patients with high severity scores who started treatment when older than 50 years, suffered a fracture within 10 years, a result consistent with that recently reported by Stiles^[25]. None of the 12 patients who started ERT during childhood or adolescence suffered a fracture during the 10 years from baseline. Additionally, regardless of the age at which treatment was begun, no patient with a low severity score had a fracture 10 years from baseline (although one patient in the > 50 age group did have a fracture event 16 years after starting ERT).

ROC analysis showed that both the baseline DS3 score and the GRAF score moderately predicted fracture risk (AUC ~0.76, 0.66 respectively) [Figure 8, Supplementary Tables 2 and 3]. For patients starting treatment between ages 18-49 years, the combined logistic model incorporating both DS3 and GRAF improved discrimination to an AUC of 0.77 (95%CI, 0.73-0.86), representing a mean ΔAUC of +0.04-0.09, which was statistically significant by bootstrap comparison (P < 0.05). Among patients aged ≥ 50 years, the DS3 and GRAF models yielded AUCs of 0.74 (95%CI, 0.66-0.82) and 0.62 (95%CI, 0.61-0.79), respectively. The combined model showed modest improvement (AUC = 0.76, 95%CI, 0.70-0.85), though the gain over single predictors was smaller (P ≈ 0.10). No reliable estimates could be derived for patients who began therapy before age 18, due to the absence of fracture events. Thus, both DS3 and GRAF individually predict fracture risk with moderate accuracy, and in combo, modestly enhance predictive performance, particularly in adults treated before age 50.

Figure 8. Receiver Operating Characteristic (ROC) curves for DS3, GRAF, and Combined (DS3 + GRAF) models predicting fracture risk during enzyme replacement therapy. Panels (A-C) depict ROC curves for each age-at-treatment group: A: < 18 years; B: 18-49 years; and C: ≥ 50 years. Each panel compares the predictive performance of the DS3 score (continuous), GRAF score (ordinal 1-3), and the combined logistic model incorporating both measures. The area under the curve (AUC) is shown for each model. The combined model modestly improved discrimination relative to single predictors, particularly in adults who began therapy before age 50, as indicated by the bootstrap-derived ΔAUC and P-values. Dashed diagonal lines represent the reference line for random classification. GRAF: Gaucher risk assessment for fracture; DS3: disease severity scoring system.

ROC/AUC shows how well a score discriminates between patients who do and do not fracture. It does not directly account for follow-up time, censoring, or whether a score predicts fracture independently of age. To address this, time-to-event analyses were performed using Cox proportional hazards regression with age included as a continuous variable. Here, baseline DS3 remained independently associated with fracture risk (HR, 1.21 per 1-point increase; 95%CI, 1.08-1.35; P = .0006). Age was also independently associated with fracture risk in this model (HR, 1.06 per year; 95%CI, 1.03-1.08; P < 0.001). In a separate model including age and GRAF, GRAF was significantly associated with fracture risk (HR, 2.33 per unit increase; 95%CI, 1.31-4.16; P = 0.004), whereas age was no longer statistically significant. This finding is consistent with the fact that age is incorporated directly into the GRAF score, and age and GRAF were strongly correlated (r = 0.84), indicating substantial overlap between these predictors. Direct comparison of hazard ratios between DS3 and GRAF should be interpreted cautiously because the scores are measured on different numerical scales and capture different underlying constructs. DS3 primarily reflects Gaucher disease severity across multiple organ systems, whereas GRAF incorporates age and history of splenectomy and bone crises to estimate fracture risk. Thus, differences in HR magnitude do not necessarily indicate superior predictive performance.

Discrimination of the Cox models was assessed using Harrell’s C-index. The age-adjusted DS3 model demonstrated a C-index of 0.80, compared with 0.78 for the age-adjusted GRAF model. Inclusion of both scores together resulted in only minimal improvement in discrimination (C-index 0.80), suggesting limited incremental predictive value of combining the two measures. Thus, pre-treatment baseline DS3 is a modestly stronger predictor of fracture risk in patients on ERT than GRAF score (KM survival and ERT), GD1 severity itself contributes to skeletal fragility beyond the effect of aging alone (COX), and inclusion of both scores together resulted in only minimal improvement in discrimination (ROC, C-index comparison).

DISCUSSION

CONCLUSIONS

Our findings suggest that DS3 captures clinically relevant aspects of Gaucher disease severity that contribute to fracture risk independently of age, supporting its potential utility as a disease-specific marker of skeletal risk in patients with Gaucher disease. In clinical practice, DS3 may be preferred for baseline risk stratification, while GRAF may retain value in settings where full DS3 assessment is not feasible. Combined use of both measures may slightly improve predictive accuracy in patients who start treatment between ages 18-50 years. However, in clinical practice, either DS3 or GRAF in conjunction with the other assessments mentioned above is likely sufficient for assessment of fracture risk. These tools may assist clinicians in identifying higher-risk patients at the time of treatment initiation although DS3, whose total score reflects severity of hematologic, splenic, hepatic, and pulmonary involvement is better suited for that purpose. However, as many other fracture risk variables are not included either in DS3 or GRAF, both should be interpreted as complementary rather than stand-alone predictive tools to assist experienced physicians in making clinical judgements and recommendations. In that regard, scientific evidence confirms that, at least for the present, based both on accumulated personal experience and appreciation of contextual skills, the ability of expert specialist physicians to assess probabilities and identify patients at risk for complications and progressive disease still trumps that of algorithm-educated, large language model artificial intelligence^[57].

Our findings support reassessment of DXA-centric monitoring strategies and highlight the need for comprehensive skeletal risk evaluation in GD1 and further study of the applicability of published guidelines for pharmacologic and non-pharmacologic treatments to prevent fractures specifically to patients with GD1^[6,58-61]. Thus, even as NBS and earlier diagnosis becomes more prevalent, it is crucial that children, teenagers and adults in whom treatment is not started at the time of diagnosis because of insufficient severity of current manifestations or failure to conform with treatment guidelines, have regular and lifetime serial expert-driven evaluations so that should disease progression occur, the risk of eventual complications, skeletal or otherwise can be minimized. Additional study topics should include patient reported and professional assessments about post-fracture functional recovery and health-related quality of life, efficacy of oral eliglustat versus ERT for prevention of GD1-associated fragility fractures, and critical evaluations of adjunctive interventions such as pharmacotherapies for osteoporosis, vitamin D optimization, and targeted weight-bearing exercises.